「Pythonではじめるベイズ機械学習入門」を読む　～第3章回帰モデル②(ポアソン回帰モデル)～

はじめに

Pythonによるベイズモデルの実装をきちんと学ぼうと思い，森賀新・木田悠歩・須山敦志著「Pythonではじめるベイズ機械学習入門」を読むことにした。

本記事は，第3章「回帰モデル」のうち一般化線形モデル(ポアソン回帰モデル)に関する読書メモである。

本書の紹介ページ

www.kspub.co.jp

3.3 一般化線形モデル : ポアソン回帰モデル

線形回帰モデルは，誤差の確率分布がガウス分布だったが，誤差を他の確率分布に拡張した回帰モデルを一般化線形モデル(generalized linear model)と呼ぶ。
本節では，誤差の分布にポアソン分布を仮定し，自然数の観測データを扱うことができるポアソン回帰モデルを紹介している。

3.3.1. モデル概要

一般化線形モデルでは，確率分布と逆リンク関数(inverse link function)を組合わせて，入力 $x$ と出力 $y$ の関係を表す。具体例としては，入力 $x$ の線形モデル $w_1x+w_2$ を逆リンク関数 $f^{-1}$ で変換して，出力 $y$ の確率分布 $p(y|\theta)$ のパラメータ $\theta$ とする。

$\begin{align} y_n &\sim p(y_n | \theta_n) \\ \\ \theta_n &= f^{-1}(w_1 x + w_2) \\ \\ f(\theta)_n &= w_1 x + w_2 \\ \end{align}$

逆リンク関数をもちいることで，線形モデル $w_1x+w_2$ を $\theta$ の定義域に合わせて変換することができる。

一般化線形モデルの例として，ポアソン回帰モデル (Poisson Regression Model)は以下のように定式化される。

$\begin{align} y_n | \lambda_n &\sim \mathrm{Poi}(y_n | \lambda_n) = \frac{ \lambda_n^{y_n} \exp(-\lambda_n) }{ y_n! } \\ \\ \lambda_n &= \exp (w_1 x_n + w_2) \in \mathbb{R}^{+} \\ \end{align}$

リンク関数に $\log$ ，すなわち逆リンク関数に $\exp$ を使うことで，ポアソン分布のパラメータ $\lambda_n$ を正の値に限定している。
またポアソン分布を用いることで，出力 $y$ を0以上の整数に限定している。

3.3.2. 実装

サンプルコードを動かしながら，挙動を確認した。
github.com

データの準備

回帰係数の真の値は， $w_1 = 0.8, w_2 = 1.2$ である。サンプルコードでは，20個の入力 $x$ について，

$\log$ リンク関数を用いて，パラメータ $\lambda \geq 0$ を得る
パラメータ $\lambda$ を用いてポアソン分布に従う乱数 $y$ を得る

という手順で，人工データを作成する。具体的な値は以下のようになっている。

出力

x_data : 
[-0.49786797  1.32194696 -2.99931375 -1.18600456 -2.11946466 -2.44596843
 -1.88243873 -0.92663564 -0.61939515  0.2329004  -0.48483291  1.111317
 -1.7732865   2.26870462 -2.83567444  1.02280506 -0.49617119  0.35213897
 -2.15767837 -1.81139107]
y_data : 
[ 6  8  0  1  3  1  0  2  0  1  1  6  0 23  0 10  5  7  0  1]

可視化すると以下のようになる。

PyMCライブラリのインポート

本書ではpymc3で説明しているが，2024年12月現在，動作確認はpymcで行なった。

サンプルコード

#import pymc3 as pm
import pymc as pm

モデルの定義

ポアソン回帰のモデル式を構築する。係数パラメータ $w_1, w_2$ にはそれぞれ事前分布を設定する必要があるが，事前分布はガウス分布を設定している。

$\begin{align} y_n | \lambda_n &\sim \mathrm{Poi}(y_n | \lambda_n) = \frac{ \lambda_n^{y_n} \exp(-\lambda_n) }{ y_n! } \\ \\ \lambda_n &= \exp (w_1 x_n + w_2) \\ \\ w_1 &\sim \mathcal{N}(0, \sigma_{w1}^2) \\ \\ w_2 &\sim \mathcal{N}(0, \sigma_{w2}^2) \\ \end{align}$

PyMCにおいて，モデルの定義は以下のように記述する。このコードにおいて， $w_1, w_2$ の事前分布は $\mathcal{N}(0, 1^2)$ としている。

サンプルコード

# モデルの定義
with pm.Model() as model:
    # 説明変数
    x = pm.Data("x", x_data)
    # 推論対象のパラメータ事前分布
    w1 = pm.Normal('w1', mu=0.0, sigma=1.0)
    w2 = pm.Normal('w2', mu=0.0, sigma=1.0)
    # 尤度関数
    y = pm.Poisson('y', mu =pm.math.exp(w1*x+w2), observed=y_data)

MCMCによる事後分布からのサンプリング

PyMCを用いて，MCMC(具体的にはNUTS)による事後分布からのサンプリングを行なった。チェーン数は3としている。チェーン数3回分のサンプルの分布とトレースプロットは下図の通りである。

3回とも定常分布に落ち着いているように見える。

次に，事後分布の可視化を行なった。

係数パラメータの推定値 $\hat{w}_1, \hat{w}_2$ のサンプル平均はそれぞれ $\hat{w}_1,=0.79, \hat{w}_2=1.2$ であった。真値である $w_1 = 0.8, w_2 = 1.2$ は，事後分布の90%信用区間の間には収まっていた。

予測分布の算出

予測分布を求めるサンプルコードは，線形回帰のときと同様に修正する必要があった。

コード(修正前)

# 検証用データ
x_new = np.linspace(-3, 3, 10)

with model:
  # 検証用データをモデルへ入力
  pm.set_data({"x": x_new})
  # 予測分布からサンプリング
  pred = pm.sample_posterior_predictive(trace, samples=1000, random_seed=1)

y_pred_samples = pred['y']

コード(修正後)

# 検証用データ
x_new = np.linspace(-3, 3, 10)

# 新しいモデルを作成
with pm.Model() as pred_model:
    # 説明変数
    x = pm.Data("x", x_new) # 検証用データで初期化

    # 事前分布 (学習済みモデルから)
    # Convert xarray DataArray to NumPy array using .values
    w1 = pm.Normal('w1', mu=trace.posterior['w1'].mean().values, sigma=trace.posterior['w1'].std().values)
    w2 = pm.Normal('w2', mu=trace.posterior['w2'].mean().values, sigma=trace.posterior['w2'].std().values)

    # 尤度関数 (観測値なし)
    y = pm.Poisson('y', mu =pm.math.exp(w1*x+w2))

    # 予測分布からサンプリング
    pred = pm.sample_posterior_predictive(trace, var_names=['y'], random_seed=1)

y_pred_samples = pred.posterior_predictive['y']

予測分布からのサンプルについても，コードを修正する。y_pred_samples の形式が(3, 3000, 10)となっている。(チェーン数，MCMCのサンプル数，予測対象の入力サンプル数)という意味なので，最初のチェーンからサンプルを持ってくるように， y_pred_samples[i,:]を y_pred_samples[0,:][i]のように修正する。

コード

# 予測分布からのサンプルを一部プロット
for i in range(0, 1000, 10):
    plt.plot(x_new, y_pred_samples[0,:][i], alpha=0.1, \
             zorder=i+1, color='green')
# 予測分布からのサンプルの平均値をプロット
plt.plot(x_new, y_pred_samples[0,:].mean(axis=0), alpha=1.0, \
         label='poisson regression', zorder=i+1, color='red')
# データ点をプロット
plt.scatter(x_data, y_data, marker='.', label='sample_data', zorder=i+2)

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.tight_layout();