pystanでポアソン回帰モデル

2020.09.02

実践Data Scienceシリーズ RとStanではじめる ベイズ統計モデリングによるデータ分析入門 (KS情報科学専門書)

「RとStanで始めるベイズ統計モデリングによるデータ分析入門」「実践編第8章 ポアソン回帰モデル」を対象に，公開されているR，Stanのコードをpython，pystanのコードへと書き直した一例です。Stanの代わりにpystanを利用しています。
本ページでは公開されていない書籍の内容については一切触れません。理論や詳しい説明は書籍を参照してください。
なお，こちらで紹介しているコードには誤りが含まれる可能性があります。内容やコードについてお気づきの点等ございましたら，ご指摘いただけると幸いです。
（本ページにて紹介しているコードはgithubにて公開しています。）

DataAnalOji.hatena.sample/python_samples/stan/3-8-ポアソン回帰モデル.ipynb at master · Data-Anal-Ojisan/DataAnalOji.hatena.sample
samples for my own blog. Contribute to Data-Anal-Ojisan/DataAnalOji.hatena.sample development by creating an account on GitHub.
github.com

分析の準備

パッケージの読み込み

import pystan
import numpy as np
import pandas as pd

%matplotlib inline
import matplotlib.pyplot as plt
plt.style.use('ggplot')
plt.rcParams['font.family'] = 'Meiryo'
import seaborn as sns

データの読み込みと図示

分析対象のデータ

fish_num_climate = pd.read_csv('3-8-1-fish-num-1.csv')
fish_num_climate.head(n=3)

データの要約

fish_num_climate.describe(include='all')

図示

# グラフ描画領域の作成
plt.figure(figsize=(10, 5))
# 散布図の描画
sns.scatterplot(x='temperature', y='fish_num', data=fish_num_climate, hue='weather')
# タイトルの変更
plt.title('釣獲尾数と気温・天気の関係')
# グラフの描画
plt.show()

brmsによる分散分析モデルの推定

pythonにbrmsの代用可能なパッケージが知っている範囲内で存在しないため，brmsを利用したコードは省略しています。

補足：brmsを用いない実装の方法

参考：デザイン行列の作成

# 応答変数の削除とダミー変数化処理
design_mat = pd.get_dummies(fish_num_climate.drop(['fish_num'], axis=1),
                   drop_first=True)

# 列名を書籍準拠に
design_mat = design_mat[['weather_sunny', 'temperature']]

# (Intercept)列追加
design_mat.insert(0, '(Intercept)', 1)

display(design_mat)

参考：データの作成

data_list_1 = dict(N=len(fish_num_climate),
                   fish_num=fish_num_climate['fish_num'],
                   temp=fish_num_climate['temperature'],
                   sunny=design_mat['weather_sunny'])

display(data_list_1)
{'N': 100,
 'fish_num': 0     0
 1     2
 2     5
 3     1
 4     3
      ..
 95    0
 96    0
 97    0
 98    0
 99    0
 Name: fish_num, Length: 100, dtype: int64,
 'temp': 0      5.5
 1     21.1
 2     17.2
 3      5.0
 4     28.3
       ... 
 95     5.9
 96    12.8
 97     2.8
 98     3.5
 99    13.2
 Name: temperature, Length: 100, dtype: float64,
 'sunny': 0     0
 1     0
 2     0
 3     0
 4     0
      ..
 95    1
 96    1
 97    1
 98    1
 99    1
 Name: weather_sunny, Length: 100, dtype: uint8}

参考：自分で変換処理を入れる

# stanコードの記述
stan_code = '''
data {
  int N;                   // サンプルサイズ
  int fish_num[N];         // 釣獲尾数
  vector[N] temp;          // 気温データ
  vector[N] sunny;         // 晴れダミー変数
}

parameters {
  real Intercept;      // 切片
  real b_temp;         // 係数(気温)
  real b_sunny;        // 係数(晴れの影響)
}

model {
  vector[N] lambda = exp(Intercept + b_temp*temp + b_sunny*sunny);
  fish_num ~ poisson(lambda);
}

'''

# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)

# サンプリング
glm_pois_stan_exp = stan_model.sampling(data=data_list_1, seed=2, n_jobs=1)

結果の表示

print(glm_pois_stan_exp.stansummary(probs=[0.025, 0.5, 0.975]))
Inference for Stan model: anon_model_10b93b39bc827a53de74eb5431c131b8.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

            mean se_mean     sd   2.5%    50%  97.5%  n_eff   Rhat
Intercept  -0.79  6.0e-3   0.24   -1.3  -0.79  -0.35   1628    1.0
b_temp      0.08  2.4e-4   0.01   0.06   0.08    0.1   1682    1.0
b_sunny    -0.59  3.7e-3   0.17  -0.93  -0.59  -0.26   2089    1.0
lp__      -37.73    0.03   1.22 -40.95 -37.41 -36.32   1412    1.0

Samples were drawn using NUTS at Tue Sep  1 19:30:40 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

参考：poisson_log関数を使用

# stanコードの記述
stan_code = '''
data {
  int N;                   // サンプルサイズ
  int fish_num[N];         // 釣獲尾数
  vector[N] temp;          // 気温データ
  vector[N] sunny;         // 晴れダミー変数
}

parameters {
  real Intercept;      // 切片
  real b_temp;         // 係数(気温)
  real b_sunny;        // 係数(晴れの影響)
}

model {
  vector[N] lambda = Intercept + b_temp*temp + b_sunny*sunny;
  fish_num ~ poisson_log(lambda);
}

generated quantities {
    vector[N] lambda_sunny;
    vector[N] lambda_cloudy;
    vector[N] fish_num_sunny;
    vector[N] fish_num_cloudy;
    
    for(i in 1:N){
        lambda_sunny[i] = Intercept +  b_temp*temp[i] + b_sunny*1;
        lambda_cloudy[i] = Intercept +  b_temp*temp[i] + b_sunny*0;
        fish_num_sunny[i] = poisson_log_rng(lambda_sunny[i]);
        fish_num_cloudy[i] = poisson_log_rng(lambda_cloudy[i]);
    }
}
'''

# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)

# サンプリング
glm_pois_stan = stan_model.sampling(data=data_list_1, n_jobs=1)

参考：結果の表示

Inference for Stan model: anon_model_6d78a2071bf51cbe1e42c5c01664ea40.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

                        mean se_mean     sd    2.5%     50%   97.5%  n_eff   Rhat
Intercept              -0.79  6.1e-3   0.24   -1.28   -0.79   -0.32   1553    1.0
b_temp                  0.08  2.4e-4   0.01    0.06    0.08     0.1   1755    1.0
b_sunny                -0.59  3.9e-3   0.17   -0.93   -0.59   -0.26   1939    1.0
lambda_sunny[1]        -0.93  4.4e-3   0.21   -1.35   -0.92   -0.54   2248    1.0
・・・省略・・・
lambda_sunny[100]      -0.29  3.1e-3   0.16   -0.61   -0.29    0.01   2596    1.0
lambda_cloudy[1]       -0.34  4.9e-3   0.19   -0.72   -0.34    0.03   1557    1.0
・・・省略・・・
lambda_cloudy[100]       0.3  3.2e-3   0.13    0.04     0.3    0.55   1690    1.0
fish_num_sunny[1]       0.41    0.01   0.65     0.0     0.0     2.0   3641    1.0
・・・省略・・・
fish_num_sunny[100]     0.74    0.01   0.86     0.0     1.0     3.0   3788    1.0
fish_num_cloudy[1]      0.73    0.01   0.86     0.0     1.0     3.0   4072    1.0
・・・省略・・・
fish_num_cloudy[100]    1.35    0.02   1.17     0.0     1.0     4.0   3936    1.0
lp__                  -37.78    0.03   1.25  -40.99  -37.46  -36.31   1352    1.0

Samples were drawn using NUTS at Tue Sep  1 19:31:19 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

ポアソン回帰の回帰曲線

95%ベイズ信頼区間付きのグラフ

MCMCサンプルの抽出

mcmc_sample = glm_pois_stan.extract()

事後予測分布を天気別にデータフレームへまとめる

result_lambda_cloudy = pd.DataFrame(np.zeros([len(fish_num_climate), 3])) 
result_lambda_sunny = pd.DataFrame(np.zeros([len(fish_num_climate), 3])) 

# 2.5・50・97.5パーセンタイルの値をデータフレームに代入
for i in range(len(fish_num_climate)):
    result_lambda_cloudy.iloc[i,:] = np.percentile(mcmc_sample['lambda_cloudy'][:,i], q=[2.5, 50, 97.5])
    result_lambda_sunny.iloc[i,:] = np.percentile(mcmc_sample['lambda_sunny'][:,i], q=[2.5, 50, 97.5])
    
# データフレームの列名を変更
result_lambda_cloudy.columns = ["2.5%", "50%", "97.5%"]
result_lambda_sunny.columns = ["2.5%", "50%", "97.5%"]

# 温度情報をデータフレームに追加
result_lambda_cloudy['temperature'] = fish_num_climate['temperature']
result_lambda_sunny['temperature'] = fish_num_climate['temperature']

# 温度情報で並び順をソート
result_lambda_cloudy = result_lambda_cloudy.sort_values('temperature')
result_lambda_sunny = result_lambda_sunny.sort_values('temperature')

回帰直線の描画

# グラフ描画領域の作成
plt.figure(figsize=(10, 5))

# 散布図の描画
sns.scatterplot(x='temperature',
                y='fish_num',
                data=fish_num_climate,
                hue='weather',
                hue_order=['cloudy', 'sunny'],
                palette=['red', 'blue'])

# 曇りの日のビールの売り上げの平均値を描画
plt.plot(result_lambda_cloudy['temperature'],
         np.exp(result_lambda_cloudy["50%"]),
         color='red')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_lambda_cloudy['temperature'],
                 y1=np.exp(result_lambda_cloudy["2.5%"]),
                 y2=np.exp(result_lambda_cloudy["97.5%"]),
                 color='red',
                 alpha=0.3)

# 晴れの日のビールの売り上げの平均値を描画
plt.plot(result_lambda_sunny['temperature'],
         np.exp(result_lambda_sunny["50%"]),
         color='blue')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_lambda_sunny['temperature'],
                 y1=np.exp(result_lambda_sunny["2.5%"]),
                 y2=np.exp(result_lambda_sunny["97.5%"]),
                 color='blue',
                 alpha=0.3)

# グラフタイトル追加
plt.title('ポアソン回帰曲線：信用区間')

# グラフの描画
plt.show()

99%ベイズ予測区間付きのグラフ

事後予測分布を天気別にデータフレームへまとめる

result_fish_num_cloudy = pd.DataFrame(np.zeros([len(fish_num_climate), 3])) 
result_fish_num_sunny = pd.DataFrame(np.zeros([len(fish_num_climate), 3])) 

# 0.5・50・99.5パーセンタイルの値をデータフレームに代入
for i in range(len(fish_num_climate)):
    result_fish_num_cloudy.iloc[i,:] = np.percentile(mcmc_sample['fish_num_cloudy'][:,i], q=[0.5, 50, 99.5])
    result_fish_num_sunny.iloc[i,:] = np.percentile(mcmc_sample['fish_num_sunny'][:,i], q=[0.5, 50, 99.5])
    
# データフレームの列名を変更
result_fish_num_cloudy.columns = ["0.5%", "50%", "99.5%"]
result_fish_num_sunny.columns = ["0.5%", "50%", "99.5%"]

# 温度情報をデータフレームに追加
result_fish_num_cloudy['temperature'] = fish_num_climate['temperature']
result_fish_num_sunny['temperature'] = fish_num_climate['temperature']

# 温度情報で並び順をソート
result_fish_num_cloudy = result_fish_num_cloudy.sort_values('temperature')
result_fish_num_sunny = result_fish_num_sunny.sort_values('temperature')

回帰直線の描画

# グラフ描画領域の作成
plt.figure(figsize=(10, 5))

# 散布図の描画
sns.scatterplot(x='temperature',
                y='fish_num',
                data=fish_num_climate,
                hue='weather',
                hue_order=['cloudy', 'sunny'],
                palette=['red', 'blue'])

# 曇りの日のビールの売り上げの平均値を描画
plt.plot(result_fish_num_cloudy['temperature'],
         result_fish_num_cloudy["50%"],
         color='red')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_fish_num_cloudy['temperature'],
                 y1=result_fish_num_cloudy["0.5%"],
                 y2=result_fish_num_cloudy["99.5%"],
                 color='red',
                 alpha=0.3)

# 晴れの日のビールの売り上げの平均値を描画
plt.plot(result_fish_num_sunny['temperature'],
         result_fish_num_sunny["50%"],
         color='blue')
# 95%ベイズ信頼区間を描画
plt.fill_between(x=result_fish_num_sunny['temperature'],
                 y1=result_fish_num_sunny["0.5%"],
                 y2=result_fish_num_sunny["99.5%"],
                 color='blue',
                 alpha=0.3)

# グラフタイトル追加
plt.title('ポアソン回帰曲線：予測区間')

# グラフの描画
plt.show()

補足：デザイン行列を使ったモデルの推定

参考：Stanに渡すデータ

data_list_2 = dict(N=len(fish_num_climate),
                   K=design_mat.shape[1],
                   Y=fish_num_climate['fish_num'],
                   X=design_mat)

display(data_list_2)
{'N': 100,
 'K': 3,
 'Y': 0     0
 1     2
 2     5
 3     1
 4     3
      ..
 95    0
 96    0
 97    0
 98    0
 99    0
 Name: fish_num, Length: 100, dtype: int64,
 'X':     (Intercept)  weather_sunny  temperature
 0             1              0          5.5
 1             1              0         21.1
 2             1              0         17.2
 3             1              0          5.0
 4             1              0         28.3
 ..          ...            ...          ...
 95            1              1          5.9
 96            1              1         12.8
 97            1              1          2.8
 98            1              1          3.5
 99            1              1         13.2
 
 [100 rows x 3 columns]}

参考：MCMCの実行

# stanコードの記述
stan_code = '''
data { 
  int N;                // サンプルサイズ
  int K;                // デザイン行列の列数(説明変数の数＋１)
  int Y[N];             // 応答変数(整数型)
  matrix[N, K] X;       // デザイン行列 
} 

parameters { 
  vector[K] b;          // 切片を含む係数ベクトル
} 

model { 
  vector[N] lambda = X * b;
  Y ~ poisson_log(lambda);
} 

'''

# モデルのコンパイル
stan_model = pystan.StanModel(model_code=stan_code)

# サンプリング
glm_pois_stan_design_mat = stan_model.sampling(data=data_list_2, seed=1, n_jobs=1)

参考：結果の表示

Inference for Stan model: anon_model_8fb9bb9e4fdb448c421e95a5e990a395.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

       mean se_mean     sd   2.5%    50%  97.5%  n_eff   Rhat
b[1]  -0.81  6.4e-3   0.24   -1.3   -0.8  -0.36   1426    1.0
b[2]  -0.59  3.5e-3   0.17  -0.93  -0.59  -0.27   2284    1.0
b[3]   0.08  2.7e-4   0.01   0.06   0.08    0.1   1433    1.0
lp__ -37.69    0.03   1.22 -40.94 -37.38 -36.31   1291    1.0

Samples were drawn using NUTS at Tue Sep  1 19:32:01 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

pystanの詳細については公式ページを参照してください。

PyStan — pystan 3.9.1 documentation
pystan.readthedocs.io

コメント