cmdstanpyの環境構築と、サンプルコード（信用区間を用いた、２群の比率の差に関する推論）
pystanほどではないが、cmdstanpyでもライブラリ間のバージョン調整が必要だったので、メモ。
トレースプロットを描画しないなら環境構築でてこずらないが、
トレースプロットを描画するならarvizが必要になり、これによって、環境構築で不具合が出てくるようになった。
stanを使うときには、venvやdockerとか使って専用環境でやったほうがいいかも。
pythonは、3.9.16
各ライブラリは、以下のバージョン。
arviz 0.12.0
cmdstanpy 0.9.68
matplotlib 3.8.4
numpy 1.23.2
scipy 1.8.1
ーー
信用区間を用いた、２群の比率の差に関する推論（２群の比率の差の分布の95%信用区間がゼロをまたぐかどうか。跨いだら差がない。跨がなければ差がありその符号でどちらが勝ちかわかる。）
のサンプルコードを以下に示す。
import cmdstanpy import numpy as np import matplotlib.pyplot as plt import arviz as az # サンプルデータ: 2群の成功数と試行数 n1, y1 = 100, 60 # 群1: 100回の試行で60回成功 n2, y2 = 120, 7 # 群2: 120回の試行で7回成功 # Stanモデルの保存 stan_code = """ data { int<lower=0> n1; // 群1のサンプルサイズ int<lower=0> y1; // 群1の成功数 int<lower=0> n2; // 群2のサンプルサイズ int<lower=0> y2; // 群2の成功数 } parameters { real<lower=0, upper=1> theta1; // 群1の成功率 real<lower=0, upper=1> theta2; // 群2の成功率 } model { theta1 ~ beta(1, 1); // 群1の成功率にベータ(1, 1)の事前分布 theta2 ~ beta(1, 1); // 群2の成功率にベータ(1, 1)の事前分布 y1 ~ binomial(n1, theta1); // 群1の成功数 y2 ~ binomial(n2, theta2); // 群2の成功数 } generated quantities { real diff; // 群1と群2の成功率の差 diff = theta1 - theta2; } """ with open('two_group_proportion.stan', 'w') as f: f.write(stan_code) # CmdStanModelオブジェクトを作成 model = cmdstanpy.CmdStanModel(stan_file='two_group_proportion.stan') # データをStanに渡す data = { 'n1': n1, 'y1': y1, 'n2': n2, 'y2': y2 } # サンプリング実行 fit = model.sample(data=data, chains=4, iter_sampling=2000, iter_warmup=1000) # 結果の表示 print(fit.summary()) # 生成された成功率の差 (diff) のサンプルを取得 diff_samples = fit.stan_variable('diff') # 95%ベイズ信用区間の計算 lower_bound = np.percentile(diff_samples, 2.5) upper_bound = np.percentile(diff_samples, 97.5) print(f'95% ベイズ信用区間: [{lower_bound:.3f}, {upper_bound:.3f}]') # 群1と群2の成功率の差の分布をプロット plt.hist(diff_samples, bins=30, density=True, alpha=0.7) plt.axvline(lower_bound, color='red', linestyle='--', label=f'2.5%: {lower_bound:.3f}') plt.axvline(upper_bound, color='blue', linestyle='--', label=f'97.5%: {upper_bound:.3f}') plt.title('Difference in Proportions (theta1 - theta2)') plt.xlabel('Difference') plt.ylabel('Density') plt.legend() plt.show() # arvizを使ってトレースプロットを描画する az_data = az.from_cmdstanpy(fit) # cmdstanpyの結果をArviZ形式に変換 # diff（成功率の差）のトレースプロットを表示 az.plot_trace(az_data, var_names=["diff"]) plt.show()
