StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第６章「練習問題」
第６章「統計モデリングの視点から確率分布の紹介」
書籍の著者 松浦健太郎 先生
この記事は、テキスト第６章「統計モデリングの視点から確率分布の紹介」「練習問題」の PyMC5写経 を取り扱います。
推論結果の妥当性については自信がありません。
はじめに
StanとRでベイズ統計モデリングの紹介
この記事は書籍「StanとRでベイズ統計モデリング」（共立出版、「テキスト」と呼びます）のベイズモデルを用いて、PyMC Ver.5で「実験的」に写経する翻訳的ドキュメンタリーです。
テキストは、2016年10月に発売され、ベイズモデリングのモデル式とプログラミングに関する丁寧な解説とモデリングの改善ポイントを網羅するチュートリアル「実践解説書」です。もちろん素晴らしいです！
「アヒル本」の愛称で多くのベイジアンに愛されてきた書籍です！
テキストに従ってStanとRで実践する予定でしたが、RのStan環境を整えることができませんでした（泣）
そこでこのシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。
引用表記
この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「StanとRでベイズ統計モデリング」初版第13刷、著者 松浦健太郎、共立出版
記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！
PyMC環境の準備
Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。
6章 練習問題
インポート
### インポート
# 数値・確率計算
import numpy as np
import scipy.stats as stats
# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')
練習問題(1)
Python の 確率統計ライブラリ scipy.stats を利用して乱数を生成します。
ベルヌーイ分布に従う乱数を生成します。
### 【離散】ベルヌーイ分布に従う乱数の生成
# パラメータ設定
p = 0.3 # 確率：0≦p≦1
size = 10 # 乱数生成数
# 乱数生成
stats.bernoulli.rvs(p=p, size=size)【実行結果】
カテゴリカル分布に従う乱数を生成します。
### 【離散】カテゴリカル分布に従う乱数の生成
# パラメータ設定
theta = [0.1, 0.2, 0.25, 0.35, 0.1] # K個、合計1
K = list(range(1, len(theta)+1)) # 1~Kの連番
size = 10 # 乱数生成数
# カテゴリカル分布器の作成
categorical = stats.rv_discrete(values=(K, theta))
# 乱数生成
categorical.rvs(size=size)【実行結果】
練習問題(2)
Python の 確率統計ライブラリ scipy.stats を利用して乱数を生成します。
確率分布の順番は第６章の掲載順です。
■ 一様分布
### 【連続】一様分布に従う乱数の生成
# パラメータ設定
a = 0 # a
b = 1 # b
size = 10 # 乱数生成数
# 乱数生成
stats.uniform.rvs(loc=a, scale=b-a, size=size)【実行結果】
■ 二項分布
### 【離散】二項分布に従う乱数の生成
# パラメータ設定
n = 8 # 試行回数n≧0の整数
p = 0.4 # 成功確率0≦p≦1
size = 10 # 乱数生成数
# 乱数生成
stats.binom.rvs(n=n, p=p, size=size)【実行結果】
■ ベータ分布
### 【連続】ベータ分布に従う乱数の生成
# パラメータ設定
a = 2 # a > 0
b = 1 # b > 0
size = 10 # 乱数生成数
# 乱数生成
stats.beta.rvs(a=a, b=b, size=size)【実行結果】
■ 多項分布
### 【離散】多項分布に従う乱数の生成
# パラメータ設定
n = 5 # 試行回数n≧0
p = [0.1, 0.6, 0.3] # 成功確率pの配列：合計1
size = 10 # 乱数生成数
# 乱数生成
stats.multinomial.rvs(n=n, p=p, size=size)【実行結果】
■ ディリクレ分布
### 【連続】ディリクレ分布に従う乱数の生成
# パラメータ設定
alpha = [1, 2, 3] # alpha > 0
size = 10 # 乱数生成数
# 乱数生成
stats.dirichlet.rvs(alpha=alpha, size=size)【実行結果】
■ 指数分布
### 【連続】指数分布に従う乱数の生成 pdf = λ * exp(-λ * x)
# パラメータ設定
lam = 2 # 平均λ > 0
size = 10 # 乱数生成数
# 乱数生成
stats.expon.rvs(scale=1/lam, size=size)【実行結果】
■ ポアソン分布
### 【離散】ポアソン分布に従う乱数の生成
# パラメータ設定
lam = 2 # 平均λ > 0
size = 10 # 乱数生成数
# 乱数生成
stats.poisson.rvs(mu=1/lam, size=size)【実行結果】
■ ガンマ分布
### 【連続】ガンマ分布に従う乱数の生成
# パラメータ設定
a = 1 # a > 0
b = 2 # b > 0
size = 10 # 乱数生成数
# 乱数生成
stats.gamma.rvs(a=a, scale=1/b, size=size)【実行結果】
■ 正規分布
### 【連続】正規分布に従う乱数の生成
# パラメータ設定
mu = 2 # 平均μ
sigma = 1.5 # 標準偏差σ > 0
size = 10 # 乱数生成数
# 乱数生成
stats.norm.rvs(loc=mu, scale=sigma, size=size)【実行結果】
■ 対数正規分布
### 【連続】対数正規分布に従う乱数の生成
# パラメータ設定
mu = 2 # 平均μ
sigma = 1.5 # 標準偏差σ > 0
size = 10 # 乱数生成数
# 乱数生成
stats.lognorm.rvs(s=sigma, scale=np.exp(mu), size=size)【実行結果】
■ 多変量正規分布
### 【連続】多変量正規分布に従う乱数の生成
# パラメータ設定：２変量正規分布の例
mean = [1, 2] # 平均μの配列
sigmas = [1, 1.5] # 標準偏差σの配列：σ > 0
rho = 0.4 # 相関係数：-1≦ρ≦1
size = 10 # 乱数生成数
# 分散共分散行列の作成
cov = [[sigmas[0]**2, rho*sigmas[0]*sigmas[1]],
[rho*sigmas[0]*sigmas[1], sigmas[1]**2]]
# 乱数生成
stats.multivariate_normal.rvs(mean=mean, cov=cov, size=size)【実行結果】
■ コーシー分布
### 【連続】コーシー分布に従う乱数の生成
# パラメータ設定
mu = 2 # 平均μ
sigma = 1.5 # 標準偏差σ > 0
size = 10 # 乱数生成数
# 乱数生成
stats.cauchy.rvs(loc=mu, scale=sigma, size=size)【実行結果】
■ Studentのt分布
### 【連続】Studentのt分布に従う乱数の生成
# パラメータ設定
nu = 4 # 自由度ν > 0
mu = 2 # 平均μ
sigma = 1.5 # 標準偏差σ > 0
size = 10 # 乱数生成数
# 乱数生成
stats.t.rvs(df=nu, loc=mu, scale=sigma, size=size)【実行結果】
■ ラプラス分布
### 【連続】ラプラス分布に従う乱数の生成
# パラメータ設定
mu = 2 # 平均μ
sigma = 1.5 # 標準偏差σ > 0
size = 10 # 乱数生成数
# 乱数生成
stats.laplace.rvs(loc=mu, scale=sigma, size=size)【実行結果】
練習問題(3)
### ２つの正規分布に従う確率変数の引き算
## 正規分布乱数の生成
# 平均μ1=0、標準偏差σ1=20、サイズ2000の正規分布乱数
y1 = stats.norm.rvs(loc=50, scale=20, size=2000, random_state=1)
# 平均μ2=20、標準偏差σ2=15、サイズ2000の正規分布乱数
y2 = stats.norm.rvs(loc=20, scale=15, size=2000, random_state=1234)
## 確率変数y=y1-y2の算出
y = y1 - y2
## 平均mu1-mu2、標準偏差√(σ1²+σ2²)の正規分布の確率密度関数zの算出
x = np.linspace(min(y), max(y), 101)
z = stats.norm.pdf(x, loc=50-20, scale=np.sqrt(20**2 + 15**2))
## 描画処理
ax = plt.subplot()
# yのヒストグラム・KDE曲線の描画
sns.histplot(y, kde=True, stat='density', ec='white', label='$y$', ax=ax)
# zの折れ線グラフムの描画
ax.plot(x, z, color='tab:red', ls='--', label='$z$')
# 修飾
ax.set(title=r'$z$ ~ $Normal\ (50-20, \sqrt{20^2+15^2})$')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles=handles[::-1], labels=labels[::-1]);【実行結果】
$${y_1}$$と$${y_2}$$はそれぞれ独立に生成されている模様です。
練習問題(4)
カイ二乗分布の乱数を生成します。
### 【連続】カイ二乗分分布に従う乱数の生成
# パラメータ設定
nu = 3 # 自由度ν
size = 10 # 乱数生成数
# 乱数生成
stats.chi2.rvs(df=nu,size=size)【実行結果】
第６章 練習問題は以上です。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で７つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
Python、EXCELのサンプルコードの配布もあります。
２．実験！たのしいベイズモデリング１＆２をPyMC Ver.5で
書籍「たのしいベイズモデリング」・「たのしいベイズモデリング２」の心理学研究に用いられたベイズモデルを PyMC Ver.5で描いて分析します。
この書籍をはじめ、多くのベイズモデルは Ｒ言語＋Stanで書かれています。
PyMCの可能性を探り出し、手軽にベイズモデリングを実践できるように努めます。
身近なテーマ、イメージしやすいテーマですので、ぜひぜひPyMCで動かして、一緒に楽しみましょう！
３．実験！岩波データサイエンス1のベイズモデリングをPyMC Ver.5で
書籍「実験！岩波データサイエンスvol.1」の４人のベイジアンによるベイズモデルを PyMC Ver.5で描いて分析します。
この書籍はベイズプログラミングのイロハをざっくりと学ぶことができる良書です。
楽しくPyMCモデルを動かして、ベイズと仲良しになれた気がします。
みなさんもぜひぜひPyMCで動かして、一緒に遊んで学びましょう！
４．楽しい写経 ベイズ・Python等
ベイズ、Python、その他の「書籍の写経活動」の成果をブログにします。
主にPythonへの翻訳に取り組んでいます。
写経に取り組むお仲間さんのサンプルコードになれば幸いです🍀
５．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析をPythonとPyMC Ver.5 で実践します。
この書籍には時系列分析のテーマが盛りだくさん！
時系列分析の懐の深さを実感いたしました。
大好きなPythonで楽しく時系列分析を学びます。
６．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Pythonのコラムを不定期に綴っています。
統計・データサイエンス書籍にまつわる記事が多いです。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
７．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learnとPyTorchの教科書です。
よかったらぜひ、お試しくださいませ。
最後までお読みいただきまして、ありがとうございました。