第17章「心理療法の介入効果」のベイズモデリングをPyMC Ver.5 で
この記事は、テキスト「たのしいベイズモデリング」の第17章「心理療法の介入効果」のベイズモデルを用いて、PyMC Ver.5で「実験的」に実装する様子を描いた統計ドキュメンタリーです。
この章では、患者への治療＝介入による改善効果について、構造方程式モデリングの主要モデルの１つ「潜在成長曲線モデル」を用いて、ベイズ推論します。
小規模の集団＝小規模なデータセットへの挑戦とのことです。
今回のPyMCモデルも紆余曲折がありました（遠い目）。
さあ一緒に、PyMCのベイズモデリングの旅へ出かけましょう！
テキストの紹介、引用表記、シリーズまえがき、PyMC等のバージョン情報は、このリンクの記事をご参照ください。
テキストで使用するデータは、R・Stan等のサンプルスクリプトとともに、出版社サイトからダウンロードして取得できます。
サマリー
テキストの概要
執筆者 ： 竹林由武 先生
モデル難易度： ★★★★★ （難しい）
自己評価
評点
$$
\begin{array}{c:c:c}
実装精度 & ★★★・・& そこそこ \\
結果再現度 & ★★・・・& やや悪い \\
楽しさ & ★★★★★& 楽しい！ \\
\end{array}
$$
評価ポイント
実はこの章のモデルを棄権する予定でした。
モデル式が複雑な行列で記述されていて判読できなかったことと、「blavaan」というパッケージで難解なコードが記述されていて解読できなかったことが理由で、PyMCのモデルを書ける気が全くしなかったのです。
しかし、ブログ化を目の前にして、テキストの文章（日本語）とパスダイアグラムを補助線にしてモデル式をPyMCに「とりあえず」書いてみたところ、なんと動いたのです！やったね！
めっちゃ頑張ったので褒めて下さい！！！
工夫・喜び・反省
クイズや謎解きのような状況でPyMCモデルをたまたま書けたのであって、モデルの本質理解には程遠い状況です。
「構造方程式」「潜在曲線モデル」このあたりを理解する必要があると痛感しております。
モデルの概要
テキストの調査・実験の概要
■ 分析対象データ
20名の患者に対して認知行動療法を実施して、次の３時点で評価した症状の得点データです。
① 介入開始時点
② 介入終了時点
③ 介入終了後４か月時点
また、介入開始前に測定された感情表出抑制傾向の得点を「共変量」として取り扱います。
図示します。
■ 潜在成長曲線モデル
テキストは「構造方程式モデリング」の１モデルである「潜在成長曲線モデル」で介入の効果を分析しています。
テキストには「潜在成長曲線モデルのパスダイアグラム」という、変数間の関係を示した図が掲載されています。
テキストよりパスダイアグラムを引用させていただきます。
この図が構造方程式モデリングを読み解く鍵になると思います！
（読み解けなかったです。。。（泣））
Pythonの構造方程式モデリング用パッケージを検索すると semopy が出てきます。
きちんと勉強しようかな・・・。
テキストのモデリング
■目的変数と関心のあるパラメータ
目的変数は、３時点の得点$${\boldsymbol{Y}}$$です。２次元です！
関心のあるパラメータは、おそらく、テキストの表17.2に掲載されている10個のパラメータだと思われます。
つまり、$${\alpha_I,\ \alpha_S,\ \gamma_{01},\ \gamma_{02},\ \psi_{I,S},\ \psi_I,\ \psi_S,\ \theta_{T1},\ \theta_{T2},\ \theta_{T3}}$$です。
■ 構造方程式
何も聞かないで下さい（聞かれても答えられません！）。
ひとまず、寛大な気持ちで数式を受け容れましょう。
$$
\begin{align*}
\boldsymbol{Y} &\sim \text{MultiNormal}\ (\mu,\ \boldsymbol{\Sigma}) \\
\mu &= \boldsymbol{\Lambda}(\alpha+\boldsymbol{\Gamma}\kappa) \\
\boldsymbol{\Sigma} &= \boldsymbol{\Lambda}(\boldsymbol{\Gamma \Phi \Gamma^{\top}}+\boldsymbol{\Psi})\boldsymbol{\Lambda^{\top}}+\boldsymbol{\Theta}
\end{align*}
$$
目的変数$${\boldsymbol{Y}}$$は症状の得点（観測変数）です。
形状は (ID, Time) であり、多変量正規分布に従います。
$${\mu}$$はモデル上の観測変数の平均ベクトル（？）です。
形状は (ID, Time) であり、いろいろなパラメータを含みます。
$${\boldsymbol{\Sigma}}$$はモデル上の観測変数の共分散行列です。
形状が (Time, Time)であり、いろいろなパラメータを含みます。
かなりツヨツヨ（強強）な見た目です！
【変数説明】
$${\boldsymbol{\Lambda}}$$は因子負荷行列（パス図の 1, 1, 1, 0, 1, 2）
$${\alpha}$$はパス図の$${\alpha_I, \alpha_S}$$であり、潜在因子の平均ベクトル
$${\boldsymbol{\Psi}}$$はパス図の$${\psi_I, \psi_S, \psi_{I,S}}$$を含む潜在因子の共分散行列
$${\boldsymbol{\Theta}}$$はパス図の$${\theta_{T1}, \theta_{T2}, \theta_{T3}}$$を含む観測変数の残差共分散行列
$${\boldsymbol{\Gamma}}$$はパス図の$${\gamma_{01}, \gamma_{02}}$$であり、時間固定共変量から切片や傾き因子へのパス係数行列
$${\kappa}$$は共変量の平均ベクトル（おそらくパス図の感情表出抑制傾向の得点）
$${\boldsymbol{\Phi}}$$は共変量の共分散行列（パス図には無さそう）
■分析・分析結果
分析の仕方や分析数値はテキストの記述が正確だと思いますので、テキストの読み込みをおすすめします。
なお、構造方程式モデリングの結果を解釈するスキルが無いため、PyMCの自己流モデルの推論値を利用した分析を行っておりません。スミマセン。
PyMC実装
Let's enjoy PyMC & Python !
準備・データ確認
１．インポート
### インポート
# ユーティリティ
import pickle
# 数値・確率計算
import pandas as pd
# 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')２．データの読み込み
csvファイルをpandasのデータフレームに読み込みます。
### データの読み込み
# データの読み込み
data = pd.read_csv('datLGM.csv')
data.rename(columns={'Unnamed: 0': 'ID'}, inplace=True)
# データの表示
print('data.shape: ', data.shape)
display(data.head())【実行結果】
全20行です。１行が患者１人分のデータです。
データ項目は次のとおりです。
「介入」はおそらく治療ということだと思います。
・ID：患者を識別するID
・BaseGRIDHAMD：介入開始前の症状の得点
・PostGRIDHAMD：介入終了後の症状の得点
・FuGRIDHAMD：フォローアップ時（介入終了後４か月時点）の症状の得点
・BaseSuppression：介入開始前に測定された感情表出抑制傾向の得点
３．データの外観の確認
データの要約統計量とヒストグラムでデータの外観を確認します。
まずは要約統計量から。
### 要約統計量の表示
display(data.iloc[:, 1:].describe().round(2))【実行結果】
最初の３つの項目は時系列です。
平均値を見ると介入終了後に得点が下がり、フォローアップ時に得点が上昇しているようです。
続いて時系列折れ線グラフを描画しましょう。
pandas の plot を利用しました。
### データの可視化：折れ線グラフ
# プロット
data.set_index('ID').T.iloc[:3, :].plot(lw=0.5)
# x軸目盛りの設定
plt.xticks([0, 1, 2], labels=['介入開始前','介入終了後','フォローアップ'])
# 凡例をグラフ外に表示
plt.legend(title='ID', bbox_to_anchor=(1.05, 1.2))
# y軸ラベル・グラフタイトル
plt.ylabel('得点')
plt.title('IDごとの得点の推移');【実行結果】
20名のIDの３時点の推移です。
介入開始前から介入終了時点にかけて、全員の得点は減少しています。
一方で、介入終了後、フォローアップ時点では、得点が上昇しているIDが多い感じです。
テキストの表17.1のデータの記述統計量を算出しましょう。
### データの記述統計量の表示 ※表17.1に相当
# データフレームの列名、行名の設定
col_names = {'BaseGRIDHAMD': '介入開始前', 'PostGRIDHAMD': '介入終了後',
'FuGRIDHAMD': 'フォローアップ', 'BaseSuppression': '共変量'}
idx_names = {'mean': '平均', 'std': '標準偏差'}
# 記述統計量の計算
data_stats = data.iloc[:, 1:].describe()
# データフレームにして表示
data_stats.iloc[1:3, :].rename(columns=col_names, index=idx_names).round(2)【実行結果】
pandas の describe() で確認した内容と同じです。
テキストの図17.1のバイオリンプロットを描画しましょう。
seaborn の violinplot() を利用します。
### 測定時点ごとの得点のバイオリンプロットの描画 ※図17.1に相当
sns.violinplot(data.iloc[:, 1:4])
plt.xticks([0, 1, 2], labels=['介入開始前','介入終了後','フォローアップ'])
plt.ylim(-11, 42)
plt.ylabel('症状の得点（HAMD）');【実行結果】
得点が下がって⤵上がる⤴傾向と、フォローアップ時点にバラツキが大きくなっていることが分かります。
モデル構築
テキストの「漠然事前分布」のモデルに取り組みます。
モデルの数式表現
目指したいPyMCのモデルの雰囲気を混ぜた「なんちゃって数式」表記です。
$$
\begin{align*}
\boldsymbol{\Lambda} &= \left[ \begin{matrix} 1 & 0 \\ 1 & 1 \\ 1 & 2\end{matrix}\right] \\
\alpha &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=10,\ \text{dims}=\text{coef}) \\
\boldsymbol{\Gamma} &\sim \text{Normal}\ (\text{mu}=0,\ \text{sigma}=10,\ \text{dims}=\text{coef}) \\
\\
\phi_{\sigma} &\sim \text{Gamma}\ (\text{alpha}=1,\ \text{beta}=0.5,\ \text{dims}=\text{coef})\\
\phi_{cov} &\sim \text{Beta}\ (\text{alpha}=1,\ \text{beta}=1)\\
\boldsymbol{\Phi} &= \left[ \begin{matrix} \phi_{\sigma0}^2 & \phi_{cov} \\ \phi_{cov} & \phi_{\sigma1}^2\end{matrix}\right] \\
\\
\psi_{\sigma} &\sim \text{Gamma}\ (\text{alpha}=1,\ \text{beta}=0.5,\ \text{dims}=\text{coef})\\
\psi_{cov} &\sim \text{Beta}\ (\text{alpha}=1,\ \text{beta}=1)\\
\boldsymbol{\Psi} &= \left[ \begin{matrix} \psi_{\sigma0}^2 & \psi_{cov} \\ \psi_{cov} & \psi_{\sigma1}^2\end{matrix}\right] \\
\\
\theta_{\sigma} &\sim \text{Gamma}\ (\text{alpha}=1,\ \text{beta}=0.5,\ \text{dims}=\text{time})\\
\theta_{cov} &\sim \text{Beta}\ (\text{alpha}=1,\ \text{beta}=1, \text{shape}=3)\\
\boldsymbol{\Theta} &= \left[ \begin{matrix} \theta_{\sigma0}^2 & \theta_{cov0} & \theta_{cov1}\\ \theta_{cov0} & \theta_{\sigma1}^2 & \theta_{cov2} \\ \theta_{cov1} & \theta_{cov2} & \theta_{\sigma2}^2 \end{matrix}\right] \\
\\
\mu &= \boldsymbol{\Lambda} (\alpha + \boldsymbol{\Gamma} \kappa) \\
\boldsymbol{\Sigma} &= \boldsymbol{\Lambda} (\boldsymbol{\Gamma} \boldsymbol{\Phi} \boldsymbol{\Gamma}^{\top} + \boldsymbol{\Psi}) \boldsymbol{\Lambda}^{\top} + \boldsymbol{\Theta}\\
likelihood &\sim \text{MvNormal}\ (\text{mu}=\mu,\ \text{cov}=\boldsymbol{\Sigma}) \\
\end{align*}
$$
１．モデルの定義
数式表現を実直にモデル記述します。
### モデルの定義
# モデルの定義
with pm.Model() as model_lcm:
### データ関連定義
## coordの定義
# 実験参加者のID
model_lcm.add_coord('id', values=range(1, 21), mutable=True)
# 時間 T1～T3
model_lcm.add_coord('time1', values=range(1, 4), mutable=True)
# 時間：２次元用
model_lcm.add_coord('time2', values=range(1, 4), mutable=True)
# 切片と傾き
model_lcm.add_coord('coef1', values=['切片', '傾き'], mutable=True)
# 切片と傾き：２次元用
model_lcm.add_coord('coef2', values=['切片', '傾き'], mutable=True)
## dataの定義
# 目的変数:得点 ※shape=(20, 3)
Y = pm.ConstantData('Y',
value=data[['BaseGRIDHAMD', 'PostGRIDHAMD', 'FuGRIDHAMD']].values,
dims=('id', 'time1'))
# 感情表出抑制傾向:BaseSuppression ※shape=(20)
kappa = pm.ConstantData('kappa',
value=data['BaseSuppression'].values, dims=('id'))
### 事前分布
# Normal(0, 10)：潜在変数の切片、因子負荷、回帰係数
# Gamma(1, 0.5)：残差、潜在変数の分散
# Beta(1, 1)：相関、共分散
# Λ：パラメータ行列 ※shape=(3, 2)
LAMBDA = pm.Deterministic('LAMBDA',
pt.stacklists([[1, 0], [1, 1], [1, 2]]),
dims=('time1', 'coef1'))
# α：潜在因子の平均ベクトル α_intercept, α_slope ※shape(2)
alpha = pm.Normal('alpha', mu=0, sigma=10, dims=('coef1',))
# Γ：時間固定共変量から切片や傾き因子へのパス係数行列 γ_01, γ_02 ※shape=(2)
GAMMA = pm.Normal('GAMMA', mu=0, sigma=10, dims=('coef1',))
# Φ：共変量の共分散行列 ※shape=(2, 2)
phiSigma = pm.Gamma('phiSigma', alpha=1, beta=0.5, dims=('coef1'))
phiCov = pm.Beta('phiCov', alpha=1, beta=1)
PHI = pm.Deterministic('PHI',
pt.stacklists([[pt.pow(phiSigma[0], 2), phiCov],
[phiCov, pt.pow(phiSigma[1], 2)]]),
dims=('coef1', 'coef2'))
# Ψ：潜在変数の共分散行列 ψ_i, ψ_s, ψ_i,s ※shape=(2, 2)
psiSigma = pm.Gamma('psiSigma', alpha=1, beta=0.5, dims=('coef1'))
psiCov = pm.Beta('psiCov', alpha=1, beta=1)
PSI = pm.Deterministic('PSI',
pt.stacklists([[pt.pow(psiSigma[0], 2), psiCov],
[psiCov, pt.pow(psiSigma[1], 2)]]),
dims=('coef1', 'coef2'))
# Θ：観測変数の残差共分散行列 θ_T1, θ_T2, θ_T3 ※shape=(3, 3)
thetaSigma = pm.Gamma('thetaSigma', alpha=1, beta=0.5, dims=('time1'))
thetaCov = pm.Beta('thetaCov', alpha=1, beta=1, shape=3)
THETA = pm.Deterministic('THETA',
pt.stacklists([[pt.pow(thetaSigma[0], 2), thetaCov[0], thetaCov[1]],
[thetaCov[0], pt.pow(thetaSigma[1], 2), thetaCov[2]],
[thetaCov[1], thetaCov[2], pt.pow(thetaSigma[2], 2)]]),
dims=('time1', 'time2'))
### muの計算 ※shape=(20, 3)
ETA = pm.Deterministic('ETA',
pt.reshape(GAMMA, (-1, 1)) @ (pt.reshape(kappa, (-1, 1))).T,
dims=('coef1', 'id'))
mu = pm.Deterministic('mu', (LAMBDA @ (alpha + ETA.T).T).T,
dims=('id', 'time1'))
### SIGMAの計算 ※shape=(3, 3)
SIGMA = pm.Deterministic('SIGMA',
LAMBDA @ (GAMMA @ PHI @ GAMMA.T + PSI) @ LAMBDA.T + THETA,
dims=('time1', 'time2'))
### 尤度
# 多変量正規分布 ※shape=(20, 3)
likelihood = pm.MvNormal('likelihood', mu=mu, cov=SIGMA, observed=Y,
dims=('id', 'time1'))【モデル注釈】
coordの定義
座標に名前を付けたり、その座標が取りうる値を設定できます。
今回は次の５つを設定しました。
・患者IDの座標：名前「id」、値「ID」
・時間の座標１：名前「time1」、値「1, 2, 3」
・時間の座標２：名前「time2」、値「1, 2, 3」
・切片・傾きの座標１：名前「coef1」、値「切片, 傾き」
・切片・傾きの座標２：名前「coef2」、値「切片, 傾き」dataの定義
目的変数である３つの時点の得点$${Y}$$（列が３つ）、感情表出抑制傾向得点$${kappa}$$を設定しました。パラメータの事前分布
テキストの文章に掲載された「Rのblavaanパッケージのデフォルト値」を設定しました（設定したつもり・・・）。
大文字はテキストで太字表記された行列です。演算子$${@}$$は、内積（ドット積）計算です。
$${\mu}$$と$${\boldsymbol{\Sigma}}$$の計算式はテキストの数式をにらめっこしつつ、内積（ドット積）や足し算処理で次元不一致エラーが解消されるように転置$${^{\top}}$$で調整して、何とか動くようになりました。
尤度
多変量正規分布です。
２．モデルの外観の確認
# モデルの表示
model_lcm【実行結果】
# モデルの可視化
pm.model_to_graphviz(model_lcm)【実行結果】
角丸四角の右下に記載される「次元 dims」を参考にして、model内のベクトル・行列の演算時に次元の不一致の解消をしました。
【ワンポイント】
演算時に次元が合っているかどうかを間接的に確認する方法をご紹介します。
もっと良い方法があるかもですが、一応ご紹介。
決定論的変数 mu の計算式が無事かどうかは mu.eval() で確認できます
mu.eval()【実行結果】
次の画像の様に乱数が表示されたなら、計算式の適切さに一歩近づいています。
演算内で不一致が起きているときには、エラーメッセージを発して処理が中断されます。
次元（形状）も確認できます。
mu.eval().shape【実行結果】
shape=(20, 3) という結果を表示しました。
３．事後分布からのサンプリング
乱数生成数はテキストに合わせています。
nuts_sampler='numpyro'とすることで、numpyroをNUTSサンプラーに利用できます。
処理時間はおよそ30秒でした。
# 事後分布からのサンプリング 30秒
with model_lcm:
idata_lcm = pm.sample(draws=5000, tune=1000, chains=4, target_accept=0.9,
nuts_sampler='numpyro', random_seed=123)【実行結果】（一部）
４．サンプリングデータの確認
$${\hat{R}}$$、事後分布の要約統計量、トレースプロットを確認します。
事後分布の収束確認はテキストにならって$${\hat{R} \leq 1.1}$$としています。
### r_hat>1.1の確認
rhat_idata = az.rhat(idata_lcm)
(rhat_idata > 1.1).sum()【実行結果】
$${\hat{R} > 1.1}$$のパラメータは「0」件です。
全てのパラメータが$${\hat{R} \leq 1.1}$$であることを確認できました。
ざっくり事後分布サンプリングデータの要約統計量とトレースプロットを確認します。
### 推論データの要約統計量の表示
pm.summary(idata_lcm, hdi_prob=0.95)【実行結果】
### トレースプロットの描画
pm.plot_trace(idata_lcm)
plt.tight_layout();【実行結果】
だいたい収束している感じがします。
共分散（例：phiCov）は０～１の幅広い区間を取っており、推測値として一抹の不安が残ります。
５．事後予測
推論の適否を確認する目的で、目的変数の事後予測サンプリングを実行します。
MCMCよりも実行時間がかかりました。およそ2分40秒です。
### 事後予測サンプリング 2分40秒
with model_lcm:
idata_lcm.extend(pm.sample_posterior_predictive(idata_lcm, random_seed=123))【実行結果】
事後予測プロットを描画します。
### 事後予測のプロット
fig, ax = plt.subplots(figsize=(10, 4))
pm.plot_ppc(idata_lcm, num_pp_samples=100, random_seed=123, ax=ax)
ax.set(xlabel='得点', ylabel='density', title='事後予測プロット');【実行結果】
データの観測値を示す黒線と比べて、オレンジ点線の事後予測の平均値は、近い位置にはあるものの、観測値のひしゃげた形状に沿うことはありません。
多変量の正規分布を仮定すると、このようなオレンジ点線になるのでしょうか。
自己流PyMCモデルの推論は観測値を適切に表現できていないかもですね。
６．事後分布サンプリングデータを用いて図表の作成
テキストの表17.3、図17.3、図17.4の作成コードを掲載します。
数値には「漠然事前分布」を設定したモデルで生成した数値を用いております。
まず、テキストの表17.3です。
### 推論データの要約統計量の表示 ※表17.3に相当
var_names = ['alpha', 'GAMMA', 'PSI', 'THETA']
pm.summary(idata_lcm, hdi_prob=0.95, var_names=var_names)【実行結果】
alpha はテキストの$${\alpha_I,\ \alpha_S}$$に近似する値になりました。
GAMMA（テキストの$${\gamma}$$）もやや近似、という感じです。
PSI（テキストの$${\phi}$$）やTHETA（テキストの$${\theta}$$）はかなり乖離している感じがします。
続いて、図17.3、17.4です。
事後分布サンプリングデータから描画するパラメータの値を抽出します。
### 図17.3, 17.4のための事後分布サンプリングデータの層別
# alpha['切片']
alpha_intercept = (idata_lcm.posterior.alpha.stack(sample=('chain', 'draw'))
.sel(coef1='切片'))
# alpha['傾き']
alpha_slope = (idata_lcm.posterior.alpha.stack(sample=('chain', 'draw'))
.sel(coef1='傾き'))
# PHI['切片', '切片']
psi_intercept = (idata_lcm.posterior.PSI.stack(sample=('chain', 'draw'))
.sel(coef1='切片').sel(coef2='切片'))
# PHI['傾き', '傾き']
psi_slope = (idata_lcm.posterior.PSI.stack(sample=('chain', 'draw'))
.sel(coef1='傾き').sel(coef2='傾き'))
# GAMMA['切片']
gamma_intercept = (idata_lcm.posterior.GAMMA.stack(sample=('chain', 'draw'))
.sel(coef1='切片'))
# GAMMA['傾き']
gamma_slope = (idata_lcm.posterior.GAMMA.stack(sample=('chain', 'draw'))
.sel(coef1='傾き'))【実行結果】なし
描画します。
### 図17.3, 17.4に相当するグラフの描画
plt.figure(figsize=(8, 10))
# alpha['切片']
ax = plt.subplot(321)
sns.histplot(alpha_intercept, bins=100, kde=True, edgecolor='white',
stat='density', ax=ax)
ax.set(title='切片因子の平均 $\\alpha_I$')
ax.grid(lw=0.5)
# alpha['傾き']
ax = plt.subplot(322)
sns.histplot(alpha_slope, bins=100, kde=True, edgecolor='white',
stat='density', ax=ax);
ax.set(title='傾き因子の平均 $\\alpha_S$')
ax.grid(lw=0.5)
# PHI['切片', '切片']
ax = plt.subplot(323)
sns.histplot(psi_intercept, bins=100, kde=True, edgecolor='white',
stat='density', ax=ax)
ax.set(title='切片因子の分散 $\\psi_I$', xlim=(-5, 60))
ax.grid(lw=0.5)
# PHI['傾き', '傾き']
ax = plt.subplot(324)
sns.histplot(psi_slope, bins=100, kde=True, edgecolor='white',
stat='density', ax=ax)
ax.set(title='傾き因子の分散 $\\psi_S$', xlim=(-5, 60))
ax.grid(lw=0.5)
# GAMMA['切片']
ax = plt.subplot(325)
sns.histplot(gamma_intercept, bins=100, kde=True, edgecolor='white',
stat='density', ax=ax)
ax.set(title='共変量から切片因子への回帰係数 $\\gamma_I$', xlim=(-1, 1))
ax.grid(lw=0.5)
# GAMMA['傾き']
ax = plt.subplot(326)
sns.histplot(gamma_slope, bins=100, kde=True, edgecolor='white',
stat='density', ax=ax)
ax.set(title='共変量から傾き因子への回帰係数 $\\gamma_S$', xlim=(-1, 1))
ax.grid(lw=0.5)
plt.tight_layout();【実行結果】
漠然事前分布を設定したモデルによる描画です。
この図はテキストの図の結果と相違しておりますので、ご注意下さい。
ただ、傾向（だいたいの雰囲気）は読み取れると思います。
【反省会】
テキストは上述の表17.3と図17.3～17.5を利用して、推論結果の分析を進めています。
しかし自己流PyMCモデルに基づく分析は、推論結果がいまひとつ精度に欠けるので、この記事では深入りしないことにします。
また、テキストの分析内容を読み解くことも難しいです（泣）
「パラメータの記号」と、切片因子・傾き因子への回帰係数といった「テキストの言葉や図表」（※）とが、どうしても紐つかず、頭に入ってきません（泣）
上述の図に用いたパラメータが適切なのかどうかも不安だったりします。
ああ、構造方程式モデリングをきちんと理解したいです。。。
（※）例えば、図17.3上の「切片因子」「傾き因子」について
これらの因子はどの変数記号のことだろうか、とか、横軸のスケールに合ったパラメータは表17.3の漠然事前分布の各値のどれなのだろうか、といった疑問でいっぱいです。
７．推論データ（idata）の保存
一応、推論データを再利用する場合に備えてファイルに保存しましょう。
idataをpickleで保存します。
### idataの保存 pickle
file = r'idata_lcm_ch16.pkl'
with open(file, 'wb') as f:
pickle.dump(idata_lcm, f)読み込みコードは次のとおりです。
### idataの読み込み pickle
file = r'idata_lcm_ch16.pkl'
with open(file, 'rb') as f:
idata_lcm_load = pickle.load(f)以上で第16章は終了です。
おわりに
数式恐怖症・数学不安症
数式を見ると体が痛くなる症状があるそうです。
私は肉体的な症状は出ないですが、厳つい変数文字が長文で並ぶ数式を見ると嫌な気持ちになります（汗）
特に気になってしまうのが総乗記号$${\prod}$$です。凱旋門です。
尤度の算出時に現れます。
$$
\displaystyle \prod^N_{i=1} f(x_i ； \theta)\cdots
$$
なので、尤度・最尤法・最尤推定量とお友達になれません。。。
シリーズの記事
次の記事
前の記事
目次
ブログの紹介
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」のシリーズが生まれています。
ベイズ書籍の実践記録も掲載中です。
最後までお読みいただきまして、ありがとうございました。
この記事が参加している募集
この記事が気に入ったらサポートをしてみませんか？