単回帰モデルによるベイズ分析
2026.07改訂
拙著「いまさら聞けないPythonでデータ分析」では、PyStan2を用いたStanによるベイズ分析を説明した。ここでは、CmdStanPyによるStanの利用を試みる。CmdStanPyの簡単なインストールの説明のウェブサイトも用意した。
拙著での単回帰分析では、緯度と予想気温という単回帰分析用の標準的なタイプのデータを扱ったが、アレルギー疾患者数の経時的変化のデータも分析してみた。度数データはポアッソン分布に従う（Spiegelhalter, 2019）と説明されているので、ポアッソン分布を用いた単回帰モデルによる分析を試みた。
単回帰モデルによる分析を行うデータとして、次のものを用いる。
RawData = [ ['City', 'Temperature',
'Latitude'],
['Sapporo',
15,
43.1],
['Sendai',
20,
38.3],
['Niigata',
17,
37.9],
['Kanazawa',
18,
36.6],
['Tokyo', 24,
35.7],
['Osaka', 22,
34.7],
['Fukuoka',
22,
33.6],
['Kochi', 24,
33.6],
['Kagoshima', 23,
31.6],
['Naha',
26,
26.2]
予想最高気温をhttps://tenki.jp/から、緯度をhttps://www.geocoding.jp/から読み取ったものである。
Latitudeを独立変数Xとし、Temperatureを従属変数Yとする単回帰モデル
Y=b0+b1*X+residual
の係数を求める。独立変数Xを中心化して求める場合と、XとYのZスコアを求めて回帰分析を行う場合を試みた。ファイルは、regres_files.zipにまとめた。
回帰モデルにおいて、独立変数を中心化した
を用いると、係数
は独立変数の平均
における回帰直線の値を表す（Gelman et al., 2021）。
式（１）のモデルによる分析は、独立変数を中心化したデータに単回帰モデルを適用すればよいので、StanスクリプトはリストA１のようになる。
リストA１ 単回帰モデルのStanスクリプト（regression.stan）
data{
int N;
array[N] real Y;
array[N] real X;
}
parameters{
real b0;
real b1;
real<lower=0.0001,
upper=1000.0> sgm;
}
transformed parameters{
array[N] real mu;
for (i in 1:N)
mu[i] = b0 + b1 * X[i]; // 回帰モデル
}
model{
b0 ~ normal(0.0, 1000.0);
b1 ~ normal(0.0, 1000.0);
sgm ~ uniform(0.0001,
1000.0);
Y ~ normal(mu, sgm);
}
リストA1を用いるPythonスクリプトはリストA2のように用意した。ファイルは、regres_files.zipにまとめた。
リストA2では、次のスクリプトによって独立変数の中心化を行っている。
X_mean = np.mean(X)
X_c = X - X_mean
# Centering
リストA2のスクリプトを次のコマンドで実行する。
(stan)
****/regres_files$ python regression_xcenter.py
リストA1のスクリプトが読み込まれてビルドが行われ、MCMCサンプリングが実行される。サンプリングが終了すると、トレース図が表示される（図A1）。
図A1
トレース図は、次のスクリプトで描いている。
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace_dist(inf_data,
var_names=['b0', 'b1', 'sgm']) #, figsize=(10,7))
plt.tight_layout()
plt.show()
図A1のWindowの右上角のX印ボタンをクリックしてWindowを閉じると、図A2の散布図が表示される。
図A2
散布図は次のスクリプトで描いている。
b0s = df_sample['b0']
b1s = df_sample['b1']
for s, t in zip(b0s, b1s):
plt.plot(s, t, 'o',
color='b', fillstyle='none')
plt.xlabel('b0', fontsize=14)
plt.ylabel('b1', fontsize=14)
plt.show()
中心化を行っているので、散布図は丸くなっている。
図A2のWindowを閉じると、図A3のグラフが表示される。
図A3
各都市の（緯度、予想最高気温）の散布図と、単回帰直線である。単回帰直線の係数は、パラメータb0およびb1の事後分布の中央値を点推定値としている（Gelman et al., 2021）。係数b0=21.10は、緯度の平均値35.13度における回帰直線の値であり、b1=-0.68は緯度が１度上がると予想気温が0.68℃下がることを表している。
図A3のWindowを閉じると、スクリプトの実行終了である。
リストA2 独立変数の中心化による回帰分析（regression_xcenter.py）
from cmdstanpy import CmdStanModel
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import arviz as az
RawData = [ ['City', 'Temperature',
'Latitude'],
['Sapporo',
15,
43.1],
['Sendai',
20,
38.3],
['Niigata',
17,
37.9],
['Kanazawa',
18,
36.6],
['Tokyo', 24,
35.7],
['Osaka', 22,
34.7],
['Fukuoka',
22,
33.6],
['Kochi', 24,
33.6],
['Kagoshima', 23,
31.6],
['Naha',
26,
26.2]
]
ID = []
Y = []
X = []
N = len(RawData) - 1
for i in range(N):
ID.append(RawData[i + 1][0])
Y.append(RawData[i + 1][1])
X.append(RawData[i + 1][2])
X_mean = np.mean(X)
X_c = X - X_mean
# Centering
Data = {'N': N, 'Y': Y, 'X': X_c}
model =
CmdStanModel(stan_file="regression.stan") # コンパイル
fit = model.sample(data=Data)
# サンプリング
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace_dist(inf_data,
var_names=['b0', 'b1', 'sgm']) #,
figsize=(10,7))
plt.tight_layout()
plt.show()
df_sample = fit.draws_pd()
# Pandas DataFrame型に変換
"""
(b0, b1)の散布図
"""
b0s = df_sample['b0']
b1s = df_sample['b1']
for s, t in zip(b0s, b1s):
plt.plot(s, t, 'o',
color='b', fillstyle='none')
plt.xlabel('b0', fontsize=14)
plt.ylabel('b1', fontsize=14)
plt.show()
"""
データ（X, Y）の散布図と回帰直線
"""
b0_med = np.median(b0s)
b1_med = np.median(b1s)
X_mean = np.mean(X)
min_X = np.min(X)
max_X = np.max(X)
for nm, xv, yv in zip(ID, X, Y):
plt.plot(xv, yv, 'o', c='b')
plt.text(xv+0.1, yv+0.1, nm,
rotation=45)
plt.plot([min_X, max_X],
[b0_med + b1_med * (min_X-X_mean), b0_med + b1_med * (max_X-X_mean)],
label='$y=b0+b1*(x-\\widehat{x})$')
plt.xlabel('Latitude', fontsize=14)
plt.ylabel('Temperature', fontsize=14)
plt.title(f'b0={b0_med:.2f}, b1={b1_med:.2f}' +
"
$\\widebar{x}$" + f"={X_mean:.2f}", fontsize=16)
plt.legend(fontsize=14)
plt.show()
独立変数Xと従属変数YのZスコアXzとYzについて単回帰分析を行う。
単回帰モデルを次のように置く。
ベイズ分析のためのStanスクリプトは、リストB1のように用意した。
リストB1 Zスコアの単回帰モデル（regressZZstan）
data{
int N;
array[N] real Y;
array[N] real X;
}
parameters{
real<lower=-10,
upper=10> a;
real<lower=0.0001,
upper=100.0> sgm;
}
transformed parameters {
array[N] real mu;
for (i in 1:N) {
mu[i] = a * X[i];
}
}
model{
a ~ uniform(-10, 10);
sgm ~ uniform(0.0001,
100.0);
Y ~ normal(mu, sgm);
}
リストB1のStanスクリプトを用いて単回帰分析を行うPythonスクリプトをリストB2のように用意した。ファイルは、regres_files.zipにまとめた。
リストB2では、Zスコアへの変換を次のスクリプトで行っている。
Xz = ss.zscore(X)
Yz = ss.zscore(Y)
リストB1のStanスクリプトのコンパイルとMCMCサンプリングの実行は次のスクリプトである。
Data = {'N': N, 'Y': Yz, 'X': Xz}
model =
CmdStanModel(stan_file="regressZZ.stan") # コンパイル
fit = model.sample(data=Data,
chains=4, inits=[{'sgm':0.5}]*4)
# サンプリング
リストB2のPythonスクリプトを次のコマンドで実行する。
(stan)
****/regres_files$ python regression_xystd.py
Stanスクリプトのコンパイルが行われ、MCMCサンプリングの終了後、トレース図が表示される（図B1）。
図B1
図B1のWindowを閉じると、図B2のグラフが表示される。
図B2
パラメータaの事後分布のグラフである。
図B2のWindowを閉じると、図B3のグラフが表示される。
図B3
各都市のデータ（緯度（Zスコア）、予想最高気温（Zスコア））の散布図と回帰直線である。回帰直線は、パラメータaの事後分布の中央値を点推定値として描いている。係数a＝−0.869であるので、緯度が１SD上がると予想最高気温が0.869SD下がる。
リストB2 Zスコアの単回帰分析Pythonスクリプト（regression_xystd.py）
from cmdstanpy import CmdStanModel
import numpy as np
import scipy.stats as ss
import matplotlib.pyplot as plt
import arviz as az
import seaborn as sb
RawData = [ ['City', 'Temperature',
'Latitude'],
['Sapporo',
15,
43.1],
['Sendai',
20,
38.3],
['Niigata',
17,
37.9],
['Kanazawa',
18,
36.6],
['Tokyo', 24,
35.7],
['Osaka', 22,
34.7],
['Fukuoka',
22,
33.6],
['Kochi', 24,
33.6],
['Kagoshima', 23,
31.6],
['Naha',
26,
26.2]
]
ID = []
Y = []
X = []
N = len(RawData) - 1
for i in range(N):
ID.append(RawData[i + 1][0])
Y.append(RawData[i + 1][1])
X.append(RawData[i + 1][2])
Xz = ss.zscore(X)
Yz = ss.zscore(Y)
Data = {'N': N, 'Y': Yz, 'X': Xz}
model =
CmdStanModel(stan_file="regressZZ.stan") # コンパイル
fit = model.sample(data=Data,
chains=4, inits=[{'sgm':0.5}]*4)
# サンプリング
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
# Arviz(az)用のデータに変換
# トレースプロット
az.plot_trace_dist(inf_data)
plt.tight_layout()
plt.show()
df_sample = fit.draws_pd()
# Pandas DataFrame型に変換
sb.kdeplot(df_sample['a'])
plt.xlabel('a', fontsize=16)
plt.title(f'Med.={np.median(df_sample["a"]):.3f}',
fontsize=16)
plt.show()
"""
データ（X, Y）の散布図と回帰直線
"""
a_med = np.median(df_sample['a'])
min_X = np.min(Xz)
max_X = np.max(Xz)
for nm, xv, yv in zip(ID, Xz, Yz):
plt.plot(xv, yv, 'o', c='b')
plt.text(xv, yv, nm,
rotation=45, fontsize=14)
plt.plot([min_X, max_X],[a_med*min_X,
a_med*max_X], label=f"y={a_med:.3f}*x")
plt.xlabel('Latitude(z-score)', fontsize=14)
plt.ylabel('Temperature(z-score)',
fontsize=14)
plt.legend(fontsize=16)
plt.show()
「アレルギー疾患の推計患者数の年次推移」のグラフ（森田・足立（2025）「アレルギーの科学」、図1‐2）があったので、厚生労働省からダウンロードしたPDFファイルhttps:/www.mhlw.go.jp/content/10900000/001545663.pdf）から、このウェブサイトのプログラムを改編利用して読み取った値を基に図P1に示すExcelファイルを作成した。
図P1 単位：1000人
アレルギー症ごとにポアッソン分布に基づく単回帰モデルを次式のように設定した。
引数
はアレルギー症の種類１：アトピー性皮膚炎、２：結膜炎、３：アレルギー性鼻炎、４：喘息を表し、引数
は時期１：平成８年、・・・、７：平成26年を表す。
式（１）はポアッソン分布を表すものであるが、生起度数データであるのでポアッソン分布を設定してみた（Spiegelhalter, 2019）。
式（２）は、平成８年の値をa[i]、平成26年の値をb[i]とする単回帰直線を表す。
上のモデルに基づくStanスクリプトを、リストP1のように作成した。ファイルは、poiregfiles.zipにまとめた。
リストP1 ポアッソン単回帰モデル（poireg.stan）
data {
array[4,7] int f;
}
parameters {
array[4] real<lower=0,
upper=1> va;
array[4] real<lower=0,
upper=1> vb;
}
transformed parameters {
array[4] real a;
array[4] real b;
array[4, 7] real lmbd;
for (i in 1:4) {
a[i]
= 1 + 2_000_000*va[i];
b[i]
= 1 + 2_000_000*vb[i];
for
(j in 1:7) {
lmbd[i][j] = a[i] + (b[i] - a[i]) * (j - 1) / 6.0;
}
}
}
model {
for (i in 1:4) {
va[i] ~ uniform(0, 1);
vb[i] ~ uniform(0, 1);
for
(j in 1:7) {
f[i][j] ~ poisson(lmbd[i][j]);
}
}
}
リストP1のStanスクリプトを用いて図P1のデータのベイズ分析を行うPythonスクリプトをリストP2のように用意した。
リストP1のスクリプトファイルとリストP2のスクリプトファイルを同じフォルダに置き、そのフォルダをカレントディレクトリとしてリスト２のスクリプトファイルを次のコマンドで実行する。
(stan)
****/poiregfiles$ python poissonreg.py
実行が始まると、図P1のファイルが読み込まれ、データが1000人を単位としているので、次のコードで調整している。
Freqs = vdata[1:]*1000
調整後、Stanスクリプトのコンパイルと実行が行われ、MCMCサンプリングが終了するとトレース図が表示される（図P2）。
図P2
図P2のWindowを閉じると、図P3のグラフが表示される。
図P3
左のグラフは、アレルギー疾患者の年度における推定患者数を表している。グラフの左上の「1e6」は、縦軸の値に1,000,000を掛けることを意味している。
実線がアレルギー疾患者のデータであり、1点鎖線が単回帰モデルによる推定値である。ポアッソン分布の場合、平均値はパラメータ値
に等しく、パラメータの点推定値として事後分布の中央値をとっている（Gelman et al., 2021）。単回帰直線は、右上がりのものと右下がりのものがある。
右側のグラフは、差b[i]-a[i]の事後分布である。事後分布が原点の右側、正の領域にあるものと、原点の左側、負の領域にあるものに分かれている。アトピー性皮膚炎とアトピー性鼻炎は事後分布が右側、正の領域にあるが、これは左のグラフにおいて回帰直線の傾きが正であること、すなわち年度が進むと患者数が増加していることを表している。
これに対して、結膜炎と喘息の事後分布は、原点の左側、負の領域にあるので、左側の単回帰直線は右下がり、すなわち年度が進むと減少していることを表している。
図P3のWindowを閉じると、スクリプトの実行終了である。
なお、現在スクリプトの実行においてmatplotlibはバージョン3.10がよいようである。matplotlib 3.11では、エラーが出た（2026.07.03）。次のコマンド
conda install
matplotlib=3.10
を実行すれば、matplotlib 3.10の最新版に置き換えられる。
リストP2 ポアッソン単回帰モデルのベイズ分析Pythonスクリプト（poissonreg.py）
from cmdstanpy import CmdStanModel
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import japanize_matplotlib
import seaborn as sb
import arviz as az
data_pd = pd.read_excel('Allerg.xlsx')
print(data_pd)
allergies = data_pd.keys()[1:]
print(allergies)
vdata = data_pd.values.T
heisei = vdata[0]
Freqs = vdata[1:]*1000
# データが千人を単位としている
print(Freqs)
sm =
CmdStanModel(stan_file='poireg.stan')
fit = sm.sample(data={'f':Freqs})
print(fit.summary())
inf_data = az.from_cmdstanpy(fit)
az.plot_trace_dist(inf_data,
var_names=['a', 'b'])
plt.tight_layout()
plt.show()
fit_pd = fit.draws_pd()
a_hat = []
b_hat = []
for i in range(4):
a_hat.append(np.median(fit_pd[f"a[{i+1}]"]))
b_hat.append(np.median(fit_pd[f"b[{i+1}]"]))
print(a_hat)
print(b_hat)
pred_lmbd = []
for i in range(4):
temp_lmbd = []
for j in range(7):
t_lmbd = a_hat[i] + (b_hat[i] - a_hat[i]) * (j - 1) / 6
temp_lmbd.append(t_lmbd)
pred_lmbd.append(temp_lmbd)
pred_lmbd = np.array(pred_lmbd)
g_color = ['g', 'orange', 'b', 'm']
plt.figure(figsize=(12,5))
plt.subplot(121)
for i in range(4):
plt.plot(heisei, Freqs[i],
label=allergies[i], color=g_color[i])
plt.plot(heisei,
pred_lmbd[i], ls='-.', color=g_color[i])
plt.plot([], ls='-.', color='aqua',
label='ポアッソン単回帰')
plt.xticks(heisei)
plt.xlabel('平成')
plt.ylabel('推計患者数')
plt.legend()
plt.title('アレルギー疾患者数年次推移')
plt.tight_layout()
plt.subplot(122)
diff = []
for i in range(4):
diff.append(fit_pd[f'b[{i+1}]'] - fit_pd[f'a[{i+1}]'])
for i in range(4):
sb.kdeplot(diff[i],
color=g_color[i], label=allergies[i])
plt.title('b[i]-a[i]の事後分布', fontsize=16)
plt.legend()
plt.tight_layout()
plt.show()
Gelman, A., Carlin, J. B., Stern, H. S., Dunson, Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014) Bayesian data analysis, third edition. CRC Press
Gelman, A., Hill, J., & Vehtari,A. (2021) Regression and other sories. Cambridge University Press.
Lambert, B. (2018) A student’s guide to Bayesian statistics. SAGE
森田英明・足立剛也（編著）（2015）アレルギーの科学．講談社
岡本安晴（2019）いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican