【改】第８章「傾いた文字は正しい文字か？鏡文字か？」のベイズモデリングをPyMC Ver.5 で
プロローグ
この記事は、以前投稿した『テキスト「たのしいベイズモデリング」の第８章「傾いた文字は正しい文字か？鏡文字か？」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。』を改造する記事です。
以前の記事はこちらです。
上記記事の執筆時点では、事後分布を収束させることができず、苦い思いをしました。
ところが投稿後、事態は一変します。
著者の先生より直々に「改善案」のご助言をいただいたのです。
この記事では先生の改善案にトライして、モデリングの改造に挑戦します。
それではテイクツー、スタートです！
リスタート
この記事は、テキスト「たのしいベイズモデリング」の第８章「傾いた文字は正しい文字か？鏡文字か？」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。
この章では、通常の文字と左右逆さまの文字（鏡文字）を０～３６０度回転させて、回答者が通常の文字か鏡文字かを判断するのにかかった反応時間を混合プロセスモデルで分析します。
今回のモデルもかなり複雑であり、自己流PyMCモデルはテキストと大きく異なる結果を出力しました（３回連続の汗）
結果はさておき、楽しくPyMCでモデリングして、ベイズ推論を満喫しましょう！
テキストの紹介、引用表記、シリーズまえがき、PyMC等のバージョン情報は、このリンクの記事をご参照ください。
テキストで使用するデータは、R・Stan等のサンプルスクリプトとともに、出版社サイトからダウンロードして取得できます。
サマリー
テキストの概要
執筆者 ： 武藤拓之 先生
モデル難易度： ★★★★★ （難しい）
収束に向けて、先生にご助言いただいたこと
データを標準化する
サンプリングするパラメータ空間の範囲をある程度小さくすると、収束効率が上がり、時間も短くなることが多い変量効果に多変量正規分布を仮定する代わりに、変量効果ごとに単変量正規分布を仮定する
ex-Gaussianのパラメータの初期値に小さい値を与える
相補誤差関数の引数が大きいと尤度が発散するため
先生、ありがとうございます！
自己評価
評点
$$
\begin{array}{c:c:c}
実装精度 & ★★★★・& ほぼ納得 \\
結果再現度 & ★★・・・& いまいち \\
楽しさ & ★★★★★& 楽しい！ \\
\end{array}
$$
評価ポイント
テキストのモデルから幾つかの変更を加えることで、なんとか収束しました。
変更内容はPyMC実装の章で。モデルを変更したことで、パラメータの事後分布統計量の見方が難しくなりました（汗）
工夫・喜び・反省
以前の（改造前の）note記事をX（旧Twitter）で投稿したところ、先生より返答をいただくことができました。
記事をお読みいただけたようで、収束のための３つの改善点を教えてくださりました。
ベイズモデリングの繋がりがネットを介して広がるのって、すごくないですか！
モデルの概要
テキストの調査・実験の概要
文字の心的回転処理のイメージ図です。
文字を見てから回答までの反応時間をコンピューターで計測しています。
■回転プロセスと非回転プロセス
図の例では文字が逆さになっていて、頭の中で回転処理をする人が多いかもしれません。
この回転処理をテキストでは「回転プロセス」と呼んでいます。
一方で、文字の傾きが小さい場合、わざわざ回転させなくても文字を判別できるでしょう。
このように回転させないで文字を判別する処理をテキストでは「非回転プロセス」と呼んでいます。
■混合プロセスモデル
２つのプロセスは反応時間にかかわります。
回転プロセスを経る場合、反応時間は長くなる可能性が高く、一方で非回転プロセスの場合、反応時間は短くなる可能性が高いのです。
「混合プロセスモデル」は、文字の向きと反応時間の対応関係に「回転プロセス」と「非回転プロセス」が混合して影響していると仮定するモデルです。
テキストに５つの仮定が記載されています。
テキストのモデリング
■目的変数と関心のあるパラメータ
目的変数$${RT}$$は反応時間（ミリ秒）です。
関心のあるパラメータは次の５つのパラメータです。
最初の３つは混合プロセスのパラメータ、残りの２つは反応時間が従うと仮定する指数－正規分布のパラメータです。
$${base}$$：基本的な処理に要する時間
$${rate}$$：心的に文字を１°回転させるのに要する時間
$${expo}$$：文字の回転の増大に伴う２つのプロセスの混合率の増加を表現
$${\sigma}$$：指数－正規分布の標準偏差パラメータ
$${\lambda}$$：指数－正規分布の平均生起回数パラメータ
■混合プロセスモデル
回転プロセスと非回転プロセスの混合モデルは整理されるとあっさりした数式になります。
$${\widehat{RT}}$$は反応時間の予測値、$${Angle}$$は文字の回転角度、添字$${i}$$は実験参加者のIDです。
この数式に至る検討過程はテキストでご確認下さい。
$$
\widehat{RT}_i = base_i + rate_i \times Angle \times \left( \cfrac{Angle}{180^{\circ}}\right)^{expo_i}
$$
■指数－正規分布モデル
反応時間は一般的には正規分布に従わず、フィッティングによく用いられる理論分布が「指数－正規分布」だそうです。
正規分布のパラメータ$${\sigma}$$と指数分布のパラメータ$${\lambda}$$によって反応時間データに含まれるノイズを表現します。
$$
RT \sim \text{ExpModNormal}\ (\ \widehat{RT}_i,\ \sigma_i,\ \lambda_i\ )
$$
■階層ベイズモデル
５つのパラメータの実験参加者の個人差が正規分布に従うと仮定しています。
個人ごとのパラメータをベクトル$${\boldsymbol{local}_i}$$で表現しています。
$$
\begin{align*}
\{ base_i, rate_i, \log expo_i, \sigma_i, \log \lambda_i \} &= \boldsymbol{local}_i \\
\boldsymbol{local}_i &\sim \text{MultiNormal}\ (\boldsymbol{\mu},\ \boldsymbol{cov}) \\
\boldsymbol{\mu} &= \{ \overline{base},\ \overline{rate},\ \log \overline{expo},\ \overline{\sigma},\ \log \overline{\lambda} \} \\
\boldsymbol{cov} &= \left[
\begin{matrix}
S_{11} &\cdots &S_{15} \\ \vdots &\ddots &\vdots \\ S_{51} &\cdots &S_{55} \\
\end{matrix}
\right]
\end{align*}
$$
$${\text{MultiNormal}}$$は多変量正規分布です。
多変量正規分布のパラメータのうち、ベクトル$${\boldsymbol{\mu}}$$は５つのパラメータの平均、行列$${\boldsymbol{cov}}$$は５つのパラメータの分散共分散行列です。
■事前分布
先行研究より$${\log \overline{expo} \cong 1.00}$$であることがわかっているとのことで、標準正規分布を弱情報事前分布としています、
その他のパラメータの事前分布については、テキストでは言及されておらず、Stanスクリプトにある程度記述されています。
$$
\log \overline{expo} \sim \text{Normal}(0, 1)
$$
■分析・分析結果
分析結果はテキストに記載の図表を利用して実施して下さい。
PyMCの自己流モデルはテキストと異なるモデルとなっております。
ご留意下さい。
PyMC実装
Let's enjoy PyMC & Python !
準備・データ確認
１．インポート
### インポート
# ユーティリティ
import pickle
# 数値・確率計算
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
plt.rcParams['font.family'] = 'Meiryo'
# ワーニング表示の抑制
import warnings
warnings.simplefilter('ignore')２．データの読み込み
csvファイルをpandasのデータフレームに読み込みます。
### データの読み込み
data_orgn = pd.read_csv('data.csv')
display(data_orgn)【実行結果】
モデリングで利用する項目は以下のとおりです。
・ID：実験参加者のID ※説明変数
・RotAngle：文字の回転角度（最短角距離：0～180°） ※説明変数
・RT：反応時間 ※目的変数
ダウンロードファイルの中に、項目説明ファイルが含まれています。
データの前処理をします。
Rスクリプトの処理と同じ内容の前処理にしました。
### データの前処理 ★Rコードに準拠
# データコピー
data = data_orgn.copy()
# 不正解試行の除外
data = data[data['Err']==0]
# 外れ値試行の除外 有効反応時間 200ms<=RT<=5000ms
data = data[(data['RT'] >= 200) & (data['RT'] <= 5000)].reset_index(drop=True)
# データの表示
display(data)【実行結果】
データ件数は11205件となりました。
３．データの外観・統計量
ひとまず要約統計量を確認します。
### 要約統計量の表示
data.describe().round(2)【実行結果】
文字タイプごとに文字の傾きと反応時間（平均値）をグラフにしてみましょう。
### 文字の傾きと平均反応時間の描画
## データの準備
# 平均反応時間
means = data.groupby(['Letter', 'Type', 'Angle'])['RT'].mean().reset_index()
# 文字のリスト
letters = means['Letter'].unique()
# タイプ（通常の文字か鏡文字か）のリスト
types = ['normal', 'mirror']
## 描画
fig, ax = plt.subplots(3, 1, figsize=(8, 8), sharey=True)
# 文字ごとに描画
for i, l in enumerate(letters):
# タイプごとに線の色を変えて折れ線グラフを描画
for t in types:
# 文字・タイプで絞り込み
tmp = means[(means['Letter']==l) & (means['Type']==t)]
# 描画
ax[i].plot(tmp['Angle'], tmp['RT'], lw=1, label=f'{t}')
# axes単位の修飾
# ax[i].set_ylim(500, 1500)
ax[i].set_xticks(np.linspace(0, 360, 7))
ax[i].set_title(f'文字:{l}')
ax[i].legend(loc='upper left')
# 全体の修飾
fig.suptitle('文字の傾きと平均反応時間')
fig.supxlabel('文字の向き（度）')
fig.supylabel('反応時間（ミリ秒）')
plt.tight_layout()
plt.show()【実行結果】
２つのグラフを作成しました。
上側が全ての反応時間、下側が500～1500秒の反応時間です。
下側の方が傾向をつかみやすいかもです。
青色の正しい文字は180°を尖った頂点にして急カーブの曲線を描いている感じがします。
オレンジ色の鏡文字は青色の線と傾向は似ていますが、大きなノイズが含まれている感じがします。
モデルの実装
階層ベイズモデルの実装です。
テキストのモデルと異なる実装を行っています。
PyMCモデリングの参考例、程度のゆるい見方をお願いします。
収束に向けた改善パターン
先生のご助言を実行したところ、尤度関数に指数-正規分布（ExGaussian）を用いる場合には収束させることができませんでした。
苦肉の策として、正規分布と指数分布の混合分布で試したところ、収束することができました！
実験した組み合わせと結果を表にまとめます。
$$
\begin{array}{c|cc}
尤度関数 & 分布 & 標準化 & 初期パラメータ & 結果\\
\hline
\text{ExGaussian} & 多変量 & 未標準化 & 未設定 &収束しない\\
& 多変量 & \boldsymbol{標準化} & 未設定 &収束しない\\
& \boldsymbol{単変量} & 未標準化 & 未設定 &収束しない\\
& 多変量 & 未標準化 & \boldsymbol{設定} & 収束しない \\
\hline
\text{混合分布} & 多変量 & 未標準化 & 未設定 & 収束しない\\
& 多変量 & \boldsymbol{標準化} & 未設定 & 収束しない\\
& \boldsymbol{単変量} & 未標準化 & 未設定 & \boldsymbol{収束した！}\\
\end{array}
$$
分布
個人差を示すパラメータが従うと仮定する正規分布の種類
テキスト通り：多変量（正規分布）、収束狙い：単変量（正規分布）標準化
データ標準化を実施した場合：標準化、実施しない場合：未標準化初期パラメータ
MCMC実行の際、パラメータ初期値に小さな値を設定した場合：設定
モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。
収束に向けた改善パターンの表で収束を示している「混合分布」「単変量（正規分布）」を用いたモデルです。
テキストとは全く異なるモデルですので、ご留意下さい。
$$
\begin{align*}
\mu_{base} &\sim \text{Uniform}\ (\text{lower}=200,\ \text{upper}=1000) \\
\mu_{rate} &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=3) \\
\mu_{\log expo} &\sim \text{TruncatedNormal}\ (\text{mu}=0,\ \text{sigma}=1, \text{lower}=-2,\ \text{upper}=2) \\
\mu_{\sigma} &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=60) \\
\mu_{\log \lambda} &\sim \text{Uniform}\ (\text{lower}=-6.5,\ \text{upper}=-4.0) \\
\sigma_{base} &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
\sigma_{rate} &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=10) \\
\sigma_{\log expo} &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=3) \\
\sigma_{\sigma} &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=100) \\
\sigma_{\log \lambda} &\sim \text{Uniform}\ (\text{lower}=0,\ \text{upper}=1.5) \\
base_i &\sim \text{Normal}\ (\text{mu}=\mu_{base},\ \text{sigma}=\sigma_{base}) \\
rate_i &\sim \text{Normal}\ (\text{mu}=\mu_{rate},\ \text{sigma}=\sigma_{rate}) \\
\log expo_i &\sim \text{Normal}\ (\text{mu}=\mu_{\log expo},\ \text{sigma}=\sigma_{\log expo}) \\
\sigma_i &\sim \text{Normal}\ (\text{mu}=\mu_{\sigma},\ \text{sigma}=\sigma_{\sigma}) \\
\log \lambda_i &\sim \text{Normal}\ (\text{mu}=\mu_{\log \lambda},\ \text{sigma}=\sigma_{\log \lambda}) \\
expo_i &= \exp (\log expo_i)\\
\lambda_i &= \exp (\log \lambda_i) \\
\widehat{RT}_i &= base_i + rate_i \times Angle \times \left( \cfrac{Angle}{180^{\circ}}\right)^{expo_i}\\
w &\sim \text{Dirichlet}\ (\text{a}=[1, 1]) \\
components &= [ \text{Normal}\ (\text{mu}=\widehat{RT}_i,\ \text{sigma}=\sigma_i), \text{Exponential}\ (\text{lam}=\lambda_i)] \\
likelihood &\sim \text{Mixture}\ (\text{w}=w,\ \text{comp\_dists}=components) \\
\end{align*}
$$
１．モデルの定義
リスト・インデックス、coord、data、パラメータの事前分布、尤度、計算値をそれぞれ定義・指定します。
id_idxは説明変数IDの値から1を差引いたインデックスです。
### モデルの定義
### リスト・インデックスの作成
# 実験参加者のIDのリストとインデックスの作成
id_list = list(set(data['ID']))
id_idx = data['ID'].values - 1
# ５つのlocalパラメータのリストの作成
params_list = ['base', 'rate', 'logExpo', 'σ', 'logλ']
### モデルの定義
with pm.Model() as model:
### データ関連定義
# coordの定義
model.add_coord('data', values=data.index, mutable=True)
model.add_coord('ID', values=id_list, mutable=True) # 実験参加者ID
model.add_coord('params', values=params_list, mutable=True) # 5つのparam
# dataの定義
RT = pm.ConstantData('RT', value=data['RT'], dims='data')
angle = pm.ConstantData('angle', value=data['RotAngle'], dims='data')
idIdx = pm.ConstantData('idIdx', value=id_idx, dims='data')
### パラメータの事前分布
## 参加者間共通の平均
muBase = pm.Uniform('muBase', lower=200, upper=1000)
muRate = pm.Uniform('muRate', lower=0, upper=3)
muLogExpo = pm.TruncatedNormal('muLogExpo', mu=0, sigma=1, lower=-2, upper=2)
muSigma = pm.Uniform('muSigma', lower=0, upper=60)
muLogLam = pm.Uniform('muLogLam', lower=-6.5, upper=-4.0)
## 参加者間共通の標準偏差
sigmaBase = pm.Uniform('sigmaBase', lower=0, upper=100)
sigmaRate = pm.Uniform('sigmaRate', lower=0, upper=10)
sigmaLogExpo = pm.Uniform('sigmaLogExpo', lower=0, upper=3)
sigmaSigma = pm.Uniform('sigmaSigma', lower=0, upper=100)
sigmaLogLam = pm.Uniform('sigmaLogLam', lower=0, upper=1.5)
## 実験参加者ごとの５つのパラメータ
base = pm.Normal('base', mu=muBase, sigma=sigmaBase, dims='ID')
rate = pm.Normal('rate', mu=muRate, sigma=sigmaRate, dims='ID')
logExpo = pm.Normal('logExpo', mu=muLogExpo, sigma=sigmaLogExpo, dims='ID')
sigma = pm.Normal('sigma', mu=muSigma, sigma=sigmaSigma, dims='ID')
logLam = pm.Normal('logLam', mu=muLogLam, sigma=sigmaLogLam, dims='ID')
# 対数をもとに戻す
expo = pm.Deterministic('expo', pt.exp(logExpo), dims='ID')
lam = pm.Deterministic('lam', pt.exp(logLam), dims='ID')
#### 尤度の定義
## 反応時間RTの予測値rtPredの算出
rtPred = pm.Deterministic(
'rtPred',
base[idIdx] + rate[idIdx] * angle * pt.pow((angle/180), expo[idIdx]),
dims='data')
## 尤度：混合分布（正規分布、指数分布） 観測値:反応時間RT
w = pm.Dirichlet('w', a=np.array([1, 1]))
components = [
pm.Normal.dist(mu=rtPred, sigma=sigma[idIdx]),
pm.Exponential.dist(lam=lam[idIdx]),
]
likelihood = pm.Mixture('likelihood', w=w, comp_dists=components,
observed=RT, dims='data')
### 計算値
muExpo = pm.Deterministic('muExpo', pt.mean(expo))
muLam = pm.Deterministic('muLam', pt.mean(lam))
muLamReciprocal = pm.Deterministic('muLamReciprocal', 1 / pt.mean(lam))【モデル注釈】
coordの定義
座標に名前を付けたり、その座標が取りうる値を設定できます。
今回は次を設定しました行の座標：名前「data」、値「行インデックス」
実験参加者のIDの座標：名前「ID」、値「ID」
パラメータの座標：名前「params」、値「param_list」（5個）
dataの定義
目的変数$${RT}$$：反応時間、説明変数$${angle}$$（RotAngle）、説明変数$${idIdx}$$（IDのインデックス）を設定しました。パラメータの事前分布
モデルの数式表現のとおりに設定しました。
テキスト・Stanスクリプトに無いパラメータ値は適当に設定しました。
実験参加者ごとの５つのパラメータは、パラメータごとに単変量正規分布を仮定しました。この部分がテキスト（多変量正規分布）と異なります。
テキストのモデルでは個人ごとのパラメータ間の関連を表現できますが、こちらのモデルでは表現できません。
尤度
正規分布と指数分布の混合分布です。
$${components}$$で正規分布と指数分布の確率 dist を設定し、$${w}$$で混合比率にディリクレ分布を定義し、pm.Mixture() で混合分布を設定しています。
テキストでは指数-正規分布を用いており、相違します。
２．モデルの外観の確認
### モデルの表示
model【実行結果】
### モデルの可視化
pm.model_to_graphviz(model)【実行結果】
巨大なグラフになりました。
テキストのモデルと相違するのは大きく次の２箇所です。
・中程の「ID(12)」で囲われた個人差の５つパラメータが、それぞれ独立した単変量正規分布です。
・左下のlilelihoodが混合分布です。
３．事後分布からのサンプリング
乱数生成数はテキストより少ない 1000draw × 4 chainです。
nuts_sampler='numpyro'とすることで、numpyroをNUTSサンプラーに利用できます。
処理時間はおよそ14分でした。
### 事後分布からのサンプリング ※NUTSサンプラーにnumpyroを使用 14分
with model:
idata = pm.sample(draws=1000, tune=1000, chains=4, target_accept=0.995,
nuts_sampler='numpyro', random_seed=999)【実行結果】省略
４．サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.01}$$としています。
### r_hat>1.1の確認
rhat_idata = az.rhat(idata)
(rhat_idata > 1.01).sum()【実行結果】
$${hat{R}}>1.01}$$のパラメータは「0」件です。
全てのパラメータが$${\hat{R} \leq 1.01}$$であることを確認できました。
パラメータ等の事後分布の要約統計量です。
パラメータ種が多いので、分割して、かつ、限定的に選んで、確認します。
### 推論データの要約統計情報の表示
var_names = ['muBase', 'muRate', 'muLogExpo', 'muSigma', 'muLogLam']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)【実行結果】
こちらは５つのパラメータの参加者間の平均です。
テキストの式(8.10)に対応するものです。
### 推論データの要約統計情報の表示
var_names = ['sigmaBase', 'sigmaRate', 'sigmaLogExpo', 'sigmaSigma', 'sigmaLogLam']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)【実行結果】
こちらは５つのパラメータの参加者間の標準偏差です。
### 推論データの要約統計情報の表示
var_names = ['w']
pm.summary(idata, hdi_prob=0.95, var_names=var_names)【実行結果】
こちらは混合分布の混合比率です。
平均値は、w[0]の正規分布が 0.864、w[1]の指数分布が 0.136 です。
トレースプロットを確認しましょう。
### トレースプロットの表示
pm.plot_trace(idata, compact=False, var_names=['muBase', 'muRate', 'muLogExpo',
'muSigma', 'muLogLam'])
plt.tight_layout();【実行結果】
### トレースプロットの表示
pm.plot_trace(idata, compact=False, var_names=['sigmaBase', 'sigmaRate',
'sigmaLogExpo', 'sigmaSigma', 'sigmaLogLam'])
plt.tight_layout();【実行結果】
### トレースプロットの表示
pm.plot_trace(idata, compact=False, var_names=['w'])
plt.tight_layout();【実行結果】
上述のトレースプロットを見る限り、分布が収束しているように感じます。
では事後予測を行って可視化し、モデルの妥当性を見てみましょう。
### 事後予測の乱数データの取得 4分20秒
with model:
idata.extend(pm.sample_posterior_predictive(idata))【実行結果】
### 事後予測プロットの描画
pm.plot_ppc(idata, num_pp_samples=4000, figsize=(12, 3));【実行結果】
黒線が観測値、青線が個別の事後予測値、オレンジ点線が事後予測平均値です。
5000ミリ秒までは事後予測が観測値に近い感じがします。
5000ミリ秒を超える観測値データはありませんので、右に裾の長い事後予測データはいまいちな部分と思われます。
５．分析～テキストにならって
テキスト表8.1の要約統計量の代替コードです。
MAPに代えてEAP（平均値）を使用しました。
### 集団レベルのパラメータの事後分布に関する要約統計量（EAPと95%HDI）
# ★テキストの表8.1に対応
# 初期値設定
val_list = ['muBase', 'muRate','muExpo', 'muSigma', 'muLamReciprocal',
'sigmaBase', 'sigmaRate','sigmaLogExpo', 'sigmaSigma', 'sigmaLogLam']
col_list = ['μ_base', 'μ_rate', 'μ_expo', 'μ_σ', 'μ_1/λ',
'σ_base', 'σ_rate', 'σ_log_expo','σ_σ', 'σ_logλ']
# 95%HDIの計算
hdi_data = pm.hdi(idata, hdi_prob=0.95)
# データフレーム化
df = pd.DataFrame()
for (val, col) in zip(val_list, col_list):
# リストの各項目のEAP（平均）・95%HDI・標準偏差をデータフレームに追加
tmp = pd.DataFrame(
{col: [idata.posterior[val].mean().data, hdi_data[val][0].data,
hdi_data[val][1].data, idata.posterior[val].std()]}
)
df = pd.concat([df, tmp], axis=1)
# データフレームの調整（列名の変更、型の変更）
df = df.T.rename(columns={0: 'EAP', 1: '2.5%',2: '97.5%', 3: 'std'}
).astype(float)
# データフレームの表示
df.round(2)【実行結果】
上５つが参加者間平均値、下５つが参加者間標準偏差（個人の差の程度）です。
テキストと違って、個人差に単変量正規分布を仮定したり、尤度関数に混合分布を用いているため、パラメータ推定値をテキストと比較することはできません。
参考値として、混合比率を掛けてみましょう。
正規分布 0.864、指数分布 0.136 を上の表の各値に一律掛けます。
なお、あくまで簡易に数値の規模感を見るための方策なので、厳密な値を求めるものではありません。
### 集団レベルのパラメータの事後分布に関する要約統計量（EAPと95%HDI）・混合比率で補正
# ★テキストの表8.1に対応
# 混合比率の平均値の取得
mix_ratio_summary = pm.summary(idata, hdi_prob=0.95, var_names='w')
mix_ratio = mix_ratio_summary['mean']
# dfの値に混合比率を乗じる
df2 = df.copy()
df2.iloc[:4] = df2.iloc[:4] * mix_ratio[0]
df2.iloc[4] = df2.iloc[4] * mix_ratio[1]
df2.iloc[5:9] = df2.iloc[5:9] * mix_ratio[0]
df2.iloc[9:] = df2.iloc[9:] * mix_ratio[1]
# データフレームの表示
df2.round(2)【実行結果】
参加者全体の平均および標準偏差です。単位はミリ秒です。
rate は１°の心的回転処理時間であり約１ミリ秒かかる感じです。
base は心的回転以外の処理時間であり約５００ミリ秒かかる感じです。
最後にテキストの表8.2を模した表を作成します。
実験参加者ごとのパラメータの平均値に混合比率を乗じた値を一覧表示します。
### 実験参加者ごとのパラメータの平均値に混合比率を乗じた値の表
# ★テキストの表8.2に相当
# 事後予測サンプリングデータから各パラメータを抽出
base_means = (idata.posterior.base.stack(sample=('chain', 'draw'))
.mean(axis=1).data)
rate_means = (idata.posterior.rate.stack(sample=('chain', 'draw'))
.mean(axis=1).data)
expo_means = (idata.posterior.expo.stack(sample=('chain', 'draw'))
.mean(axis=1).data)
sigma_means = (idata.posterior.sigma.stack(sample=('chain', 'draw'))
.mean(axis=1).data)
lam_means = (idata.posterior.lam.stack(sample=('chain', 'draw'))
.mean(axis=1).data)
# データフレーム形式で表示
pd.DataFrame({'ID': list(range(1, 13)),
'base': base_means * mix_ratio[0],
'rate': rate_means * mix_ratio[0],
'expo': expo_means * mix_ratio[0],
'σ': sigma_means * mix_ratio[0],
'1/λ': 1 / lam_means * mix_ratio[1]
}).round(2)【実行結果】
１つ前の表は集団レベルの値、こちらの表は個人レベルの値です。
ざっくりと base、rate、expo、1/λ はテキストの値が少々上ブレしている感じです。
σの値はテキストとかなり相違しており、こちらの方がバラツキが大きい感じです。
おまけ的に、テキストの図8.6のグラフを作成してみます。
実験参加者ごとの反応時間の予測値$${\widehat{RT}}$$の事後分布サンプリングデータを用いて、平均値・95%HDI区間・50%HDI区間を算出し、観測値とともに描画します。
テキストと比べてHDI区間が激狭なのが気になりますが、ひとまずコードを描きます。
なお、この描画では混合比率を考慮しておりません。
まずは描画用のデータ作成から。
### ID・Angle・平均・95%HDI・50%HDIのデータフレームを作成
# rtPredの事後分布データを抽出・dataデータフレームと結合して、data_postを作成
# rtPredとID・Angleとの関係を繋ぐ
rt_pred_sample = idata.posterior.rtPred.stack(sample=('chain', 'draw')).data
data_post = pd.concat([data, pd.DataFrame(rt_pred_sample)], axis=1)
# ID・Angle・平均・95%HDI・50%HDIのデータフレームを作成
plot_data = []
# IDごとに繰り返し処理
for id in range(1, 13):
tmp_data = data_post[data_post['ID']==id]
# ID単位で角度Angleごとに計算処理を繰り返す
for angle in sorted(tmp_data['Angle'].unique()):
# ID・Angle単位のrtPredデータを平坦化して取得
post = tmp_data[tmp_data['Angle']==angle].iloc[:, 14:].values.flatten()
# ID・Angle単位の95%HDIの算出
hdi95low, hdi95high = az.hdi(post, hdi_prob=0.95)
# ID・Angle単位の50%HDIの算出
hdi50low, hdi50high = az.hdi(post, hdi_prob=0.50)
# ID・Angle単位の平均・95%HDI・50%HDIをリストに追加
plot_data.append((id, angle, post.mean(), hdi95low, hdi95high,
hdi50low, hdi50high))
# データフレーム化
plot_df = pd.DataFrame(plot_data,
columns=['ID', 'Angle', 'mean', '2.5%', '97.5%', '25%', '75%'])
display(plot_df)【実行結果】
ID、角度Angle、$${\widehat{RT}}$$の事後平均とHDI４値をデータフレーム plot_df にまとめました。
続いて描画します。
### plot_dfをIDごとに描画 ★図8.6に相当
fig, ax = plt.subplots(4, 3, sharey=True, figsize=(10, 10))
# IDごとに描画処理を繰り返す
for i in range(1, 13):
# 処理対象IDのデータの取り出し
tmp_plot = plot_df[plot_df['ID']==i]
tmp_plot2 = data[data['ID']==i]
# 描画するaxesを特定
axes = divmod((i-1),3)
# rtPredの事後平均を描画（青色）
ax[axes].plot(tmp_plot['Angle'], tmp_plot['mean'], color='tab:blue')
# rtPredの95%HDIを描画（オレンジ色）
ax[axes].fill_between(tmp_plot['Angle'], tmp_plot['2.5%'],
tmp_plot['97.5%'], color='tab:orange',)
# rtPredの50%HDIを描画（緑色）
ax[axes].fill_between(tmp_plot['Angle'], tmp_plot['25%'],
tmp_plot['75%'], color='tab:green', alpha=0.5)
# 観測値RTの描画（赤色）
ax[axes].scatter(tmp_plot2['Angle'], tmp_plot2['RT'], color='tab:red',
s=0.1)
# 修飾
ax[axes].set(ylim=(200, 2000), title=f'ID:{i}')
plt.tight_layout()
plt.show()【実行結果】
小さな赤い点が観測値、青線が$${\widehat{RT}}$$の平均値、オレンジ領域が$${\widehat{RT}}$$の95%HDI区間、緑領域が$${\widehat{RT}}$$の50%HDI区間です。
HDI区間が狭すぎる感じがします
$${\widehat{RT}}$$はMCMCサンプリング時に取得しているので、モデルが良くなかったのかな・・・
個人レベルの標準偏差$${\sigma}$$が大きいことが影響しているのかな・・・
いまひとつ腑に落ちない結末でした。
６．推論データ（idata）の保存
推論データを再利用することはないと思いますが、ひとまずファイルに保存しましょう。
idataをpickleで保存します。
### idataの保存 pickle
file = r'idata_ch08.pkl'
with open(file, 'wb') as f:
pickle.dump(idata, f)読み込みコードは次のとおりです。
### idataの読み込み pickle
file = r'idata_ch08.pkl'
with open(file, 'rb') as f:
idata_load = pickle.load(f)以上で第８章は終了です。
おわりに
鏡文字か、鏡餅か
人間の認知プロセスに切り込む認知心理学的なテーマでした。
予めベイズモデリングを定め、実験参加者から実験データを収集して、分析するという、一連のデータ分析過程を疑似体験できる素晴らしい章です。
がしかし、PyMC自己流モデルの限界を超えてしまい、事後分布（からの乱数）が収束することなく、BADエンディングを迎える結果になってしまいました。
残念です・・・。
とは言うものの、この章のモデルに取り組んだ当初は、サンプリング時のエラーが解消できず、「お蔵入り」候補だったのです。
エラー回避ができるようになったことだけでも喜んでおこうと思います！
やったー！
と書いたのが最初の投稿のときです。
その後の推移はプロローグ等のとおりです。
現時点の状況は、指数-正規分布を用いたモデルで収束できておらず、PyMCのExGaussianクラスを使いこなせていないという課題が継続しています。
しかし、先生のアドバイスが優しく背中を押してくれて、別の分布を用いて収束まで漕ぎ着けられました。
先生、本当にありがとうございます！
テキストの分析をPyMCで追いかけることについては道半ばであるものの、ベイズモデリングを調整する体験はとても貴重なものになりました。
清々しく新年を迎えられそうです！
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
note で４つのシリーズ記事を書いています。
ぜひ覗いていってくださいね！
１．のんびり統計
統計検定２級の問題集を手がかりにして、確率・統計をざっくり掘り下げるブログです。
雑談感覚で大丈夫です。ぜひ覗いていってくださいね。
統計検定２級公式問題集CBT対応版に対応しています。
２．RとStanではじめる心理学のための時系列分析入門 を PythonとPyMC Ver.5 で
書籍「RとStanではじめる心理学のための時系列分析入門」の時系列分析トピックを PythonとPyMC Ver.5で取り組みます。
豊富なテーマ（お題）を実践することによって、PythonとPyMCの基礎体力づくりにつながる、そう信じています。
日々、Web検索に勤しみ、時系列モデルの理解、Pythonパッケージの把握、R・Stanコードの翻訳に励んでいます！
このシリーズがPython時系列分析の入門者の参考になれば幸いです🍀
３．Python機械学習プログラミング実践記
書籍「Python機械学習プログラミング PyTorch & scikit-learn編」を学んだときのさまざまな思いを記事にしました。
この書籍は、scikit-learn と PyTorch の教科書です。
よかったらぜひ、お試しくださいませ。
４．データサイエンスっぽいことを綴る
統計、データ分析、AI、機械学習、Python のコラムを不定期に綴っています。
「統計」「Python」「数学とPython」「R」のシリーズが生まれています。
ベイズ書籍の実践記録も掲載中です。
最後までお読みいただきまして、ありがとうございました。
この記事が気に入ったらサポートをしてみませんか？