ポアッソン分布に基づく回帰分析

Stanによるベイズ分析

改訂2025.10

Spiegelhalter (2019)は、EnglandとWalesにおける1998年～2016年の各年ごとの殺人事件の件数を取り上げている。Spiegelhalter(2019), p. xviに記載されているウェブサイトからダウンロードしたファイルを参考にして、回帰曲線を多項式で表し、残差をポアッソン分布に従うとする一般化非線形回帰モデル（Gelmanら（2014）、Gelmanら（2021））を適用してみた。図１に示す結果を得た。

図１

横軸に年、縦軸に件数をとり、多項式とポアッソン分布に基づく回帰曲線が描かれている。Spiegelhalter (2019)は、件数にポアッソン分布を当てはめている。ポアッソン分布は、次式で表される。

このとき、

が成り立つ。

いま、パラメータが、独立変数の多項式で表される場合について考える。

すなわち、

である。ここで、独立変数の平均が０であるようにしておくと、は独立変数が0のときのの値を表す。

上のモデルによりベイズ分析を行うStanスクリプトをリスト１のように用意した。ファイルはpoisson_regression.zipにまとめた。

リスト１のStanスクリプトを用いて、多項式（３）の次数を１から６まで変えたときのモデルのデータへの適合度をLOOおよびWAICで表すスクリプトをリスト２のように用意した。CmdStanPy対応である。CmdStanPyの準備の簡単な説明を、このウェブサイトに用意した。

リスト２のスクリプトを実行すると、図２のグラフが表示される。

図２

次数pの値ごとに、多項式（３）が描かれている。

各次数のモデルの適合度は、スクリプトとの実行終了後、端末に表示されている以下のようなメッセージを見ればよい。

p = 1

loo = 305.5710445476431

waic = 306.1465984928944

p = 2

loo = 258.9537429248399

waic = 257.9747681073625

p = 3

loo = 198.32892678074032

waic = 196.49849902593593

p = 4

loo = 192.3191932158917

waic = 191.4583530362248

p = 5

loo = 194.5435322926794

waic = 193.1524949731479

p = 6

loo = 194.3558604721531

waic = 192.76692761849955

(stan) *****/poisson_regression$

loo（PSIS-LOO-CV）、waic（WAIC）、ともにｐ＝４のときが最小である。多項式は４次式が最適であるとして、ｐ＝４のときのベイズ分析のスクリプトをリスト３のように作成した。

リスト３のスクリプトを実行すると、図３のグラフが表示される。

図３

係数パラメータの事後分布の中央値を係数とする回帰曲線（式（3））のグラフが、黒の曲線で描かれている。データ値を表す棒グラフによく合っていると思われる。回帰曲線上に描かれている両端が＋記号のライム色の線分は、パラメータが黒の曲線で表される回帰曲線の値であるポアッソン分布の2.5％値と97.5％値を端点とする線分である。2005年を除いて、棒グラフの示す値はこの［2.5％値、97.5％値］の間にある。2005年の頻度は、回帰モデルで表される全体の傾向から統計的に有意に外れていると解釈される。2005年は、ロンドン同時爆破事件のあった年である。

黒の破線および1点鎖線の曲線は、MCMCサンプルからランダムに1000組を取り出し、それぞれの組の値で回帰曲線を描いて、それらの2.5％値と97.5％値を描いたものである。1000本の曲線から1／3を選んで橙色の曲線で描いている。

橙色の曲線群、および2.5％点と97.5％点を繋いだ曲線は、回帰曲線の不確定度を表すものであり、ライム色の線分は頻度値（観測値）の不確定度を表すものであると言える。

なお、子供の交通事故について、同じようにポアッソン分布を用いた非線形回帰モデルを用いた最尤法による分析を、このウェブサイトで行っている。

リスト１　ポアッソン回帰モデルのStanスクリプト（poisson_reg_dp_loo_waic.stan）

// ポアッソン回帰モデル

// 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）

from cmdstanpy import CmdStanModel

import matplotlib.pyplot as plt

import numpy as np

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) # Centering

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}

# コンパイル

model = CmdStanModel(stan_file='poisson_reg_dp_loo_waic.stan')

# サンプリング

fit = model.sample(data={'p':p, 'n':len(years_c), 'x':years_c, 'y':counts},

chains=4, inits=[init_f()]*4)

print(fit.summary())

i_data = az.from_cmdstanpy(fit,log_likelihood='LL' ) # Arviz(az)用のデータに変換

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.draws_pd()

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])

リスト３　回帰多項式を4次式としてベイズ分析を行うPythonスクリプト（p_4.py）

from cmdstanpy import CmdStanModel

import matplotlib.pyplot as plt

import numpy as np

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の初期値

# コンパイル

model = CmdStanModel(stan_file='poisson_reg_dp_loo_waic.stan')

# サンプリング

fit = model.sample(data={'p':p, 'n':len(years_c), 'x':years_c, 'y':counts},

chains=4, inits=[init_f()]*4)

print(fit.summary())

inf_data = az.from_cmdstanpy(fit) # Arviz(az)用のデータに変換

# トレース図

az.plot_trace(inf_data, var_names=['a','b','lmbd'])

plt.tight_layout()

plt.show()

d_fm = fit.draws_pd() # Pandas DataFrame

plt.figure(figsize=(9,7))

plt.bar(years, counts)

# 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()

参考文献

Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., & Rubin, D. B. (2014). Bayesian data analysis, 3rd. Ed. CRC Press.

Gelman, A., Hill, J., & Vehtari, A. (2021). Regression and other stories. Cambridge University

岡本安晴（2009）　統計学を学ぶための数学入門［下］．培風館

岡本安晴（2019）　いまさら聞けないPythonでデータ分析．丸善出版

Spiegelhalter, D. (2019). The art of statistics: Learning from data. Pelican

Home