ポアッソン分布に基づく回帰分析
実行時エラーと対応（2024.09.03）
Ubuntu 22/WSLでの実行時エラー：
「cannot import name ’gaussian’ from ‘scipy.signal’ 」というエラーが表示された。これは、scipyをバージョン1.12.0のものに置き換えると解決した。
置き換えは、インストールと同じコマンド
conda install
scipy==1.12.0
で行った。
Ubuntu 24/WSLでの実行時エラー：
Scipyに関するエラーは、Ubuntu 22の場合と同様、scipy 1.12.0をインストールすることで対応できた。
その後、実行すると
「Could not load the Qt platform plugin “xcb”」なるエラーが表示されるようになった。これは、表示デバイスがらみのエラーのようであるので、Ubuntu 24/WSLでは、StanによるMCMCサンプリングまでとして、後の処理はWindowsで行うことにした。例を、このウェブサイトの後半に示す。
Spiegelhalter (2019)は、EnglandとWalesにおける1998年〜2016年の各年ごとの殺人事件の件数を取り上げている。Spiegelhalter(2019), p. xviに記載されているウェブサイトからダウンロードしたファイルを参考にして、リスト１に示すスクリプトを作成した。リスト１のスクリプトを実行すると、図１のグラフが表示される。
図１
横軸に年、縦軸に件数をとり、ポアッソン分布に基づく回帰曲線が描かれている。Spiegelhalter (2019)は、件数にポアッソン分布を当てはめている。ポアッソン分布は、次式で表される（岡本、2009；岡本、2019）。
このとき、
が成り立つ。
いま、パラメータ
が、独立変数
の多項式で表される場合について考える。
すなわち、
である。ここで、独立変数の平均が０であるようにしておくと、
は独立変数が0のときのの値を表す。
上のモデルによりベイズ分析を行うStanスクリプトをリスト２のように用意した。
回帰モデルにおける多項式（3）の次数を求めるために、次数を１から６まで変化させて、モデルの適合度を調べるプログラムを、リスト３のように用意した。
リスト３のスクリプトを実行すると、１次多項式から６次多項式までの各回帰多項式モデルによるベイズ分析の結果が、図２のように表示される。
図２
４次多項式以上のモデルにおいて、回帰曲線はよくデータに当てはまっているように思われる。
モデル適合度の指標としてLOO-CVとWAICがリスト３のスクリプトにより出力されるが、それらの値を表１とした。
|
表１ モデルの適合度 |
||
|
多項式 |
LOO-CV |
WAIC |
|
１次式モデル |
306.54 |
305.99 |
|
２次式モデル |
260.67 |
258.33 |
|
３次式モデル |
198.03 |
196.26 |
|
4次式モデル |
192.25 |
191.32 |
|
5次式モデル |
195.16 |
193.69 |
|
6次式モデル |
194.24 |
192.81 |
モデル適合度の値から判断すると、LOO-CVおよびWAICのいずれにおいても４次多項式モデルが最も良いと考えられる。
リスト１のプログラムは、４次多項式による分析を行うものである。実行すると下図（図１再掲）の結果が得られる。
図１（再掲）
係数パラメータの事後分布の中央値を係数とする回帰曲線（式（3））のグラフが、黒の曲線で描かれている。データ値を表す棒グラフによく合っていると思われる。回帰曲線上に描かれている両端が＋記号のライム色の線分は、パラメータが黒の曲線で表される回帰曲線の値であるポアッソン分布の2.5％値と97.5％値を端点とする線分である。2005年を除いて、棒グラフの示す値はこの［2.5％値、97.5％値］の間にある。2005年の頻度は、回帰モデルで表される全体の傾向から統計的に有意に外れていると解釈される。
黒の破線および1点鎖線の曲線は、MCMCサンプルからランダムに1000組を取り出し、それぞれの組の値で回帰曲線を描いて、それらの2.5％値と97.5％値を描いたものである。1000本の曲線から1／3を選んで橙色の曲線で描いている。
橙色の曲線群、および2.5％点と97.5％点を繋いだ曲線は、回帰曲線の不確定を表すものであり、ライム色の線分は頻度値（観測値）の不確定を表すものであると言える。
リスト１ ポアッソン分布に基づく回帰分析（p_4.py、圧縮ファイルpoisson_regression.zipに含めた）
import matplotlib.pyplot as plt
import numpy as np
import stan
import arviz as az
import scipy.stats as ss
"""
データは、Spiegelhalter(2019)のウェブサイト
https://github.com/dspiegel29/ArtofStatistics
からダウンロードしたファイル09-4-homocide-trends-x.csv
に含まれているものである。
"""
years = np.array([1997, 1998, 1999,
2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012,
2013, 2014, 2015])
counts = np.array([587, 609, 648, 676,
761, 736, 729, 749, 639,
676, 712, 617, 574, 599, 497, 508, 506, 497, 557])
years_c = years - np.mean(years) # 中心化
p=4
# 4次の多項式を用いる
#
#
MCMCサンプリング
#
def init_f():
return {'v_a':
np.mean(counts), 'v_b':[0.0]*p}
# MCMCの初期値
with
open('poisson_reg_dp_loo_waic.stan') as f_stan:
sm =
stan.build(f_stan.read(),
dict({'p':p, 'n':len(years_c), 'x':years_c, 'y':counts}))
fit = sm.sample(num_chains=4,
init=[init_f()]*4)
i_data = az.from_pystan(posterior=fit,
log_likelihood='LL', posterior_model=sm)
print(az.summary(i_data))
plt.figure(figsize=(9,7))
plt.bar(years, counts)
d_fm = fit.to_frame()
#
#
4次の多項式
#
def reg_4(x, va, b1, b2, b3, b4):
return va + b1*x + b2*(x**2)
+ b3*(x**3) + b4*(x**4)
smpl_a = d_fm['a']
smpl_b1 = d_fm['b.1']
smpl_b2 = d_fm['b.2']
smpl_b3 = d_fm['b.3']
smpl_b4 = d_fm['b.4']
#
# MCMCサンプルからランダムに1000組のサンプルを選び
# それぞれのサンプル値による回帰曲線を描く
#
n_samples = len(smpl_a)
rng = np.random.default_rng()
boot_idx = rng.choice(n_samples,
1000) # ランダムに1000個採るための添字
regs = []
for i in range(1000):
idx = boot_idx[i]
reg_v = reg_4(years_c,
smpl_a[idx], smpl_b1[idx], smpl_b2[idx],
smpl_b3[idx], smpl_b4[idx])
regs.append(reg_v)
if (i % 3) == 0:
# ランダムな回帰曲線を間引いて描画
plt.plot(years, reg_v, c='orange', ls=':', alpha=0.1)
#
# ランダムに選んだ回帰曲線の2.5%値と97.5%値を描く
#
reg025_975 = np.percentile(regs, [2.5,
97.5], axis=0)
plt.plot(years, reg025_975[0], c='k',
ls='--', label='2.5% of regs')
plt.plot(years, reg025_975[1], c='k',
ls='-.', label='97.5% of regs')
#
# パラメータの事後分布の中央値を係数値とする
# 回帰曲線を描く
#
reg_v = reg_4(years_c,
np.median(smpl_a), np.median(smpl_b1),
np.median(smpl_b2), np.median(smpl_b3),
np.median(smpl_b4))
plt.plot(years, reg_v, c='k',
label='Reg.by median\ndegree = 4')
#
# 中央値を係数値とする回帰曲線に対応するポアッソン分布の
#
2.5%値と97.5%値を端点する区間を描く
#
pp025 = []
pp975 = []
for v in reg_v:
pp025.append(ss.poisson.ppf(0.025, v))
pp975.append(ss.poisson.ppf(0.975, v))
for i, v in enumerate(years):
plt.plot([v,v], [pp025[i],
pp975[i]], lw=5,
markersize=20, marker='+', c='lime')
plt.plot([], lw=5, c='lime',
label='95% interval')
plt.title('Homicides Each Year in
England and Wales', fontsize=16)
plt.xticks(years, years,
rotation='vertical')
plt.xlabel('Year', fontsize=18)
plt.ylabel('Count of Incidents',
fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()
リスト２ ポアッソン分布による回帰分析Stanスクリプト（poisson_reg_dp_loo_waic.stan、圧縮ファイルpoisson_regression.zipに含めた）
//
//
ポアッソン回帰モデル
//
Y. Okamoto, 2024.08
//
data {
int p;
// 多項式の次数
int n;
// データ数
array[n] real x; //
独立変数
array[n] int y; // 従属変数
}
transformed data {
real mean_y;
mean_y = mean(y);
array[n] vector[p] Xs; // 独立変数のj乗の値の配列
for (i in 1:n)
for
(j in 1:p) {
Xs[i][j] = x[i] ^ j;
}
}
parameters {
//
// PyStan 3の要請に合わせたパラメータ
// transformed parametersでモデルに合うパラメータに変換
//
real<lower=0.0> v_a;
vector[p] v_b;
}
transformed parameters {
real a;
// 回帰式のｙ切片
vector[p] b;
// 回帰式の偏回帰係数
array[n] real lmbd; // ポアッソン分布のパラメタ
//
// リスケーリング
//
a = mean_y * v_a;
for (j in 1:p)
b[j]
= v_b[j] * mean_y/(10^j);
//
// ポアッソン分布のパラメータlmbdの回帰式
//
for (i in 1:n) {
lmbd[i] = a;
for
(j in 1:p) {
lmbd[i] += b[j] * Xs[i][j];
if (lmbd[i] <= 0.0)
// 非正値は強制的に変更
lmbd[i] = 0.001;
}
}
}
model {
v_a ~ exponential(0.1);
v_b ~ normal(0.0, 10.0);
for (i in 1:n)
y[i]
~ poisson(lmbd[i]); // ポアッソン回帰モデル
}
generated quantities {
array[n] real LL; // 対数尤度
for (i in 1:n)
LL[i] = poisson_lpmf(y[i] | lmbd[i]);
}
リスト３ 多項次数を決めるためのプログラム（p_1to6.py、圧縮ファイルpoisson_regression.zipに含めた）
import matplotlib.pyplot as plt
import numpy as np
import stan
import arviz as az
"""
データは、Spiegelhalter(2019)のウェブサイト
https://github.com/dspiegel29/ArtofStatistics
からダウンロードしたファイル09-4-homocide-trends-x.csv
に含まれているものである。
"""
years = np.array([1997, 1998, 1999,
2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012,
2013, 2014, 2015])
counts = np.array([587, 609, 648, 676,
761, 736, 729, 749, 639,
676, 712, 617, 574, 599, 497, 508, 506, 497, 557])
for v in zip(years, counts):
print(v)
years_c = years - np.mean(years)
years_c = np.array(years_c)
years = np.array(years)
counts = np.array(counts, dtype=int)
Loos = np.empty(6)
waics = np.empty(6)
plt.figure(figsize=(12,5))
for p in range(1, 7):
#
#
MCMCサンプリング
#
def init_f():
return {'v_a': np.mean(counts), 'v_b':[0.0]*p}
with
open('poisson_reg_dp_loo_waic.stan') as f_stan:
sm =
stan.build(f_stan.read(),
dict({'p':p, 'n':len(years_c), 'x':years_c, 'y':counts}))
fit =
sm.sample(num_chains=4, init=[init_f()]*4)
i_data =
az.from_pystan(posterior=fit, log_likelihood='LL', posterior_model=sm)
print(az.summary(i_data))
Loo = az.loo(i_data,
scale='deviance')
Loos[p-1] = Loo.elpd_loo
waic = az.waic(i_data,
scale='deviance')
waics[p-1] = waic.elpd_waic
plt.subplot(int(f'23{p}'))
plt.bar(years, counts,
label=f'p={p}')
d_fm = fit.to_frame()
med_a = np.median(d_fm['a'])
med_b = np.empty(p)
for j in range(p):
med_b[j] = np.median(d_fm[f'b.{j+1}'])
def reg_p(x):
v =
med_a;
for
i in range(p):
v += med_b[i] * x**(i+1)
return v
reg_v = reg_p(years_c)
plt.plot(years, reg_v)
plt.legend()
plt.tight_layout()
plt.show()
for p in range(1,7):
print('p =', p)
print('loo =', Loos[p-1])
print('waic =', waics[p-1])
岡本安晴（2009） 統計学を学ぶための数学入門［下］．培風館
岡本安晴（2019） いまさら聞けないPythonでデータ分析．丸善出版
Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican
Ubuntu 24/WSLでQtがらみのエラーが表示されるので（2024.09.03）、以下のようにプログラムをUbuntu 24/WSLで実行する部分（リストUW１）と、その結果を受けてPython 12/Windows 10で実行する部分（UW2）に分けて実行した。
リストUW1の実行後、リストUW2を実行すると、図UW1が表示される。
図UW1
#import matplotlib.pyplot as plt
import numpy as np
import stan
import arviz as az
import scipy.stats as ss
import pickle
"""
データは、Spiegelhalter(2019)のウェブサイト
https://github.com/dspiegel29/ArtofStatistics
からダウンロードしたファイル09-4-homocide-trends-x.csv
に含まれているものである。
"""
years = np.array([1997, 1998, 1999,
2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012,
2013, 2014, 2015])
counts = np.array([587, 609, 648, 676,
761, 736, 729, 749, 639,
676, 712, 617, 574, 599, 497, 508, 506, 497, 557])
years_c = years - np.mean(years) # 中心化
p=4
# 4次の多項式を用いる
#
#
MCMCサンプリング
#
def init_f():
return {'v_a':
np.mean(counts), 'v_b':[0.0]*p}
# MCMCの初期値
with
open('poisson_reg_dp_loo_waic.stan') as f_stan:
sm =
stan.build(f_stan.read(),
dict({'p':p, 'n':len(years_c), 'x':years_c, 'y':counts}))
fit = sm.sample(num_chains=4,
init=[init_f()]*4)
i_data = az.from_pystan(posterior=fit,
log_likelihood='LL', posterior_model=sm)
print(az.summary(i_data))
with open('i_data.pkl', 'wb') as f:
pickle.dump(i_data, f)
#plt.figure(figsize=(9,7))
#plt.bar(years, counts)
d_fm = fit.to_frame()
with open('d_fm.pkl', 'wb') as f:
pickle.dump(d_fm, f)
リストUW2 Windowsでの実行スクリプト（リストUW1の実行後に実行）
import matplotlib.pyplot as plt
import numpy as np
#import stan
import arviz as az
import scipy.stats as ss
import pickle
"""
データは、Spiegelhalter(2019)のウェブサイト
https://github.com/dspiegel29/ArtofStatistics
からダウンロードしたファイル09-4-homocide-trends-x.csv
に含まれているものである。
"""
years = np.array([1997, 1998, 1999,
2000, 2001, 2002, 2003, 2004,
2005, 2006, 2007, 2008,
2009, 2010, 2011, 2012,
2013, 2014, 2015])
counts = np.array([587, 609, 648, 676,
761, 736, 729, 749, 639,
676, 712, 617, 574, 599, 497, 508, 506, 497, 557])
years_c = years - np.mean(years) # 中心化
p=4
# 4次の多項式を用いる
plt.figure(figsize=(9,7))
plt.bar(years, counts)
with open('i_data.pkl', 'rb') as f:
i_data = pickle.load(f)
with open('d_fm.pkl', 'rb') as f:
d_fm = pickle.load(f)
#
# 4次の多項式
#
def reg_4(x, va, b1, b2, b3, b4):
return va + b1*x + b2*(x**2)
+ b3*(x**3) + b4*(x**4)
smpl_a = d_fm['a']
smpl_b1 = d_fm['b.1']
smpl_b2 = d_fm['b.2']
smpl_b3 = d_fm['b.3']
smpl_b4 = d_fm['b.4']
#
# MCMCサンプルからランダムに1000組のサンプルを選び
# それぞれのサンプル値による回帰曲線を描く
#
n_samples = len(smpl_a)
rng = np.random.default_rng()
boot_idx = rng.choice(n_samples,
1000) # ランダムに1000個採るための添字
regs = []
for i in range(1000):
idx = boot_idx[i]
reg_v = reg_4(years_c,
smpl_a[idx], smpl_b1[idx], smpl_b2[idx],
smpl_b3[idx], smpl_b4[idx])
regs.append(reg_v)
if (i % 3) == 0:
# ランダムな回帰曲線を間引いて描画
plt.plot(years, reg_v, c='orange', ls=':', alpha=0.1)
#
# ランダムに選んだ回帰曲線の2.5%値と97.5%値を描く
#
reg025_975 = np.percentile(regs, [2.5,
97.5], axis=0)
plt.plot(years, reg025_975[0], c='k',
ls='--', label='2.5% of regs')
plt.plot(years, reg025_975[1], c='k',
ls='-.', label='97.5% of regs')
#
# パラメータの事後分布の中央値を係数値とする
# 回帰曲線を描く
#
reg_v = reg_4(years_c,
np.median(smpl_a), np.median(smpl_b1),
np.median(smpl_b2), np.median(smpl_b3),
np.median(smpl_b4))
plt.plot(years, reg_v, c='k',
label='Reg.by median\ndegree = 4')
#
#
中央値を係数値とする回帰曲線に対応するポアッソン分布の
#
2.5%値と97.5%値を端点する区間を描く
#
pp025 = []
pp975 = []
for v in reg_v:
pp025.append(ss.poisson.ppf(0.025, v))
pp975.append(ss.poisson.ppf(0.975, v))
for i, v in enumerate(years):
plt.plot([v,v], [pp025[i],
pp975[i]], lw=5,
markersize=20, marker='+', c='lime')
plt.plot([], lw=5, c='lime',
label='95% interval')
plt.title('Homicides Each Year in
England and Wales', fontsize=16)
plt.xticks(years, years,
rotation='vertical')
plt.xlabel('Year', fontsize=18)
plt.ylabel('Count of Incidents',
fontsize=14)
plt.legend()
plt.tight_layout()
plt.show()