StanとRでベイズ統計モデリングをPyMC Ver.5で写経～第11章「11.2 混合正規分布」

第11章「離散値をとるパラメータを使う」

書籍の著者　松浦健太郎 先生

---

この記事は、テキスト第11章「離散値をとるパラメータを使う」の 11.2節「混合正規分布」の PyMC5写経 を取り扱います。
PyMCの混合正規分布クラス NormalMixture() を用います。

はじめに

---

StanとRでベイズ統計モデリングの紹介

この記事は書籍「StanとRでベイズ統計モデリング」（共立出版、「テキスト」と呼びます）のベイズモデルを用いて、PyMC Ver.5で「実験的」に写経する翻訳的ドキュメンタリーです。

StanとRでベイズ統計モデリング (Wonderful R 2)

テキストは、2016年10月に発売され、ベイズモデリングのモデル式とプログラミングに関する丁寧な解説とモデリングの改善ポイントを網羅するチュートリアル「実践解説書」です。もちろん素晴らしいです！
「アヒル本」の愛称で多くのベイジアンに愛されてきた書籍です！

テキストに従ってStanとRで実践する予定でしたが、RのStan環境を整えることができませんでした（泣）
そこでこのシリーズは、テキストのベイズモデルをPyMC Ver.5に書き換えて実践します。

引用表記

この記事は、出典に記載の書籍に掲載された文章及びコードを引用し、適宜、掲載文章とコードを改変して書いています。
【出典】
「StanとRでベイズ統計モデリング」初版第13刷、著者 松浦健太郎、共立出版

記事中のイラストは、「かわいいフリー素材集いらすとや」さんのイラストをお借りしています。
ありがとうございます！

かわいいフリー素材集

PyMC環境の準備

Anacondaを用いる環境構築とGoogle ColaboratoryでPyMCを動かす方法について、次の記事にまとめています。
「PyMCを動かすまでの準備」章をご覧ください。

11.2 混合正規分布

---

インポート

### インポート

# 数値・確率計算
import pandas as pd
import numpy as np

# PyMC
import pymc as pm
import pytensor.tensor as pt
import arviz as az

# 描画
import matplotlib.pyplot as plt
import seaborn as sns
plt.rcParams['font.family'] = 'Meiryo'

# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')

モデル式11.5

２つの正規分布が混合したと想定されるデータを取り扱います。

データの読み込み・確認

サンプルコードのデータを読み込みます。

### データの読み込み ◆データファイル11.3 data-mix1.txt
# Y: 能力測定スコア

data1 = pd.read_csv('./data/data-mix1.txt')
print('data1.shape: ', data1.shape)
display(data1.head())

【実行結果】

データの外観を確認します。
Y のヒストグラムとKDE曲線を描画します。
テキスト図11.2に相当します。

### ヒストグラムとKDE曲線の描画 ◆図11.2

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))
twinx1 = ax.twinx()  # KDE曲線の軸
# ヒストグラムの描画
sns.histplot(data=data1, x='Y', bins=25, ec='white', ax=ax)
# KDE曲線の描画
sns.kdeplot(data=data1, x='Y', fill=True, alpha=0.2, ax=twinx1)
# 修飾
ax.grid(lw=0.5)
twinx1.set(ylim=(0, 0.22), yticks=[], ylabel='');

【実行結果】
峰が２つあるように見えます。

PyMCのモデル定義

PyMCでモデル式11-5を実装します。周辺化消去はしません。
モデルの定義です。
ラベルスイッチング対策には muDiff を用いています。
muDiff≦0, mu1 + muDiff = mu1 とすることで、mu1≦mu2 を表現します。 

### モデルの定義 ◆モデル式11-5 model11-5.stan

# 混合正規分布の構成数の設定
n_components = 2

# モデルの定義
with pm.Model() as model1:
    
    ### データ関連定義
    ## coordの定義
    model1.add_coord('data', values=data1.index, mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data1['Y'].values, dims='data')

    ### 事前分布
    # 混合比率 weight
    a = pm.Uniform('a', lower=0, upper=1)
    weight = pt.stack([a, 1-a])
    # 混合正規分布の平均パラメータμ ※mudiffはラベルスイッチング対策
    mu1 = pm.Uniform('mu1', lower=-10, upper=10)
    muDiff = pm.Uniform('muDiff', lower=0, upper=10)
    mu2 = pm.Deterministic('mu2', mu1 + muDiff)
    mu = pt.stack([mu1, mu2])
    # 混合正規分布の標準偏差パラメータσ
    sigma = pm.Uniform('sigma', lower=0, upper=10, shape=n_components)

    ### 尤度関数 混合正規分布
    obs = pm.NormalMixture('obs', w=weight, mu=mu, sigma=sigma, observed=Y,
                           dims='data')

モデルの定義内容を見ます。

### モデルの表示
model1

【実行結果】

### モデルの可視化
pm.model_to_graphviz(model1)

【実行結果】

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 10秒
with model1:
    idata1  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.8,
                        nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

### r_hat>1.1の確認
# 設定
idata_in = idata1        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
print((az.rhat(idata_in) > threshold).sum())

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['a', 'mu1', 'mu2', 'muDiff', 'sigma']
pm.summary(idata1, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
pm.plot_trace(idata1, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

推定結果の解釈

事後分布の要約統計量を算出します。
算出関数を定義します。

### mean,sd,2.5%,25%,50%,75%,97.5%パーセンタイル点をデータフレーム化する関数の定義
def make_stats_df(y):
    probs = [2.5, 25, 50, 75, 97.5]
    columns = ['mean', 'sd'] + [str(s) + '%' for s in probs]
    quantiles = pd.DataFrame(np.percentile(y, probs, axis=0).T, index=y.columns)
    tmp_df = pd.concat([y.mean(axis=0), y.std(axis=0), quantiles], axis=1)
    tmp_df.columns=columns
    return tmp_df

要約統計量を算出します。

### 要約統計量の算出・表示
vars = ['a', 'mu1', 'mu2']
param_samples = idata1.posterior[vars].to_dataframe().reset_index(drop=True)
sigma_samples = pd.DataFrame(
            idata1.posterior.sigma.stack(sample=('chain', 'draw')).data.T,
            columns=[f'sigma[{i+1}]' for i in range(2)])
param_df = pd.concat([param_samples, sigma_samples], axis=1)
display(make_stats_df(param_df).round(3))

【実行結果】
テキスト２１０ページに掲載の各パラメータの中央値・95%ベイズ信頼区間とほぼ同じ結果を得られました。

事後分布プロットを描画します。

### パラメータの事後分布の描画
pm.plot_posterior(idata1, hdi_prob=0.95, point_estimate='median', round_to=3,
                  var_names=['a', 'mu1', 'mu2', 'sigma'],
                  grid=(2, 3), figsize=(10, 7))
plt.tight_layout();

【実行結果】

モデル式11.6

複数（数は未定）の正規分布が混合したと想定されるデータを取り扱います。

データの読み込み・確認

サンプルコードのデータを読み込みます。

### データの読み込み ◆データファイル11.4 data-mix2.txt
# Y: 能力測定スコア

data2 = pd.read_csv('./data/data-mix2.txt')
print('data2.shape: ', data2.shape)
display(data2.head())

【実行結果】

データの外観を確認します。
Y のヒストグラムとKDE曲線を描画します。
テキスト図11.3に相当します。

### ヒストグラムとKDE曲線の描画 ◆図11.3

# 描画領域の設定
fig, ax = plt.subplots(figsize=(5, 4))
twinx1 = ax.twinx()  # KDE曲線の軸
# ヒストグラムの描画
sns.histplot(data=data2, x='Y', bins=25, ec='white', ax=ax)
# KDE曲線の描画
sns.kdeplot(data=data2, x='Y', fill=True, alpha=0.2, ax=twinx1)
# 修飾
ax.grid(lw=0.5)
twinx1.set(ylim=(0, 0.16), yticks=[], ylabel='');

【実行結果】
峰が３つあるように見えます。

PyMCのモデル定義

PyMCでモデル式11-6を実装します。周辺化消去はしません。
モデルの定義です。
ラベルスイッチング対策には muDiff を用いています。
muDiff≦0, mu1 + muDiff = mu1 とすることで、mu1≦mu2 を表現します。 
混合する正規分布の数を５としており、for文を回して５つの mu を定義しています。

### モデルの定義 ◆モデル式11-6 model11-6.stan

# 混合正規分布の構成数の設定
n_components = 5  # テキストのK

# モデルの定義
with pm.Model() as model2:
    
    ### データ関連定義
    ## coordの定義
    model2.add_coord('data', values=data2.index, mutable=True)
    ## dataの定義
    # 目的変数 Y
    Y = pm.ConstantData('Y', value=data2['Y'].values, dims='data')

    ### 事前分布

    ## 混合比率 weight
    weight = pm.Dirichlet('weight', a=np.ones(n_components), shape=n_components)

    ## 混合正規分布の平均パラメータμ shape=n_components
    ## ※muDiff[k]はラベルスイッチング対策
    # mu0の標準偏差
    sigmaMu = pm.Uniform('sigmaMu', lower=0, upper=50)
    # muの差分muDiff0～3の標準偏差
    sigmaMuDiff = pm.Uniform('sigmaMuDiff', lower=0, upper=50)
    # for文で回す際のmu, muDiffの受け皿listの作成
    mu = [0] * n_components
    muDiff = [0] * (n_components - 1)
    # muの0番目の事前分布
    mu[0] = pm.Normal('mu0', mu=data2['Y'].mean(), sigma=sigmaMu)
    # muDiffの0～3番目の事前分布とmu（mu[k]+muDiff[k]）の算出
    for k in range(0, n_components-1):
        muDiff[k] = pm.HalfNormal('muDiff'+str(k), sigma=sigmaMuDiff)
        mu[k+1] = pm.Deterministic('mu'+str(k+1), mu[k] + muDiff[k])
    # muをひとまとめにする
    mu = pt.stack([mu[k] for k in range(n_components)])

    ## 混合正規分布の標準偏差パラメータσ shape=n_components
    sigma = pm.Gamma('sigma', alpha=1.5, beta=1.0, shape=n_components)

    ### 尤度関数 混合正規分布
    obs = pm.NormalMixture('obs', w=weight, mu=mu, sigma=sigma, observed=Y,
                           dims='data')

モデルの定義内容を見ます。

### モデルの表示
model2

【実行結果】

### モデルの可視化
pm.model_to_graphviz(model2)

【実行結果】
mu0 ～ mu4 の５つの mu の定義が複雑になってしまいました。

MCMCの実行と収束確認

MCMCを実行します。

### 事後分布からのサンプリング 10秒
with model2:
    idata2  = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.9,
                        nuts_sampler='numpyro', random_seed=1234)

【実行結果】省略

Pythonで事後分布からのサンプリングデータの確認を行います。
Rhatの確認から。
テキストの収束条件は「chainを3以上にして$${\hat{R}<1.1}$$のとき」です。

### r_hat>1.1の確認
# 設定
idata_in = idata2        # idata名
threshold = 1.01         # しきい値

# しきい値を超えるR_hatの個数を表示
display((az.rhat(idata_in) > threshold).sum())

【実行結果】
収束条件を満たしています。

事後統計量を表示します。

### 推論データの要約統計情報の表示
var_names = ['mu0', 'mu1', 'mu2', 'mu3', 'mu4', 'sigmaMu', 'sigmaMuDiff',
             'weight', 'sigma', 'muDiff0', 'muDiff1', 'muDiff2', 'muDiff3']
pm.summary(idata2, hdi_prob=0.95, var_names=var_names, round_to=3)

【実行結果】

トレースプロットを描画します。

### トレースプロットの表示
pm.plot_trace(idata2, compact=True, var_names=var_names)
plt.tight_layout();

【実行結果】

推定結果の解釈

事後分布の要約統計量を算出します。

### 要約統計量の算出・表示
vars1 = ['mu0', 'mu1', 'mu2', 'mu3', 'mu4', 'sigmaMu', 'sigmaMuDiff']
vars2 = ['muDiff0', 'muDiff1', 'muDiff2', 'muDiff3']
param_samples1 = idata2.posterior[vars1].to_dataframe().reset_index(drop=True)
param_samples2 = idata2.posterior[vars2].to_dataframe().reset_index(drop=True)
weight_samples = pd.DataFrame(
            idata2.posterior.weight.stack(sample=('chain', 'draw')).data.T,
            columns=[f'weight[{i+1}]' for i in range(5)])
sigma_samples = pd.DataFrame(
            idata2.posterior.sigma.stack(sample=('chain', 'draw')).data.T,
            columns=[f'sigma[{i+1}]' for i in range(5)])
param_df = pd.concat([param_samples1, weight_samples, sigma_samples,
                      param_samples2], axis=1)
display(make_stats_df(param_df).round(3))

【実行結果】
テキストに推論結果が掲載されていないため、このモデルの推論の適否は分かりません・・・。

事後分布プロットを描画します。

### パラメータの事後分布の描画
pm.plot_posterior(idata2, hdi_prob=0.95, point_estimate='median', round_to=3,
                  var_names=['mu0', 'mu1', 'mu2', 'mu3', 'mu4', 'weight',
                             'sigma', 'sigmaMu', 'sigmaMuDiff'],
                  grid=(4, 5), figsize=(12, 12))
plt.tight_layout();

【実行結果】

11.2 節は以上です。
なお、model11-6b.stanはStan特有の文法と思われるため、PyMC化は省略しました。

---

シリーズの記事

次の記事

前の記事

目次

ブログの紹介

---

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の教科書です。
よかったらぜひ、お試しくださいませ。

最後までお読みいただきまして、ありがとうございました。