🦔
cmdstanpy : 重回帰
import os
import time
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from cmdstanpy import CmdStanModel
データの読み込み
以下のデータを使用しています。
df = pd.read_csv("./data-attendance-1.txt")
print(df.head())
A Score Y
0 0 69 0.286
1 1 145 0.196
2 0 125 0.261
3 1 86 0.109
4 1 158 0.230
- Scoreは正規分布に従いそう
- A=1の人はYが小さそう
- Scoreが高いとYも高そう
sns.pairplot(df, hue="A")
<seaborn.axisgrid.PairGrid at 0x28b499ac0>
モデル式
上記を見るとAとScoreの線型結合がYに線形に相関すると考えて、線形の重回帰でモデル式を構築します。
output_dir = "./output/"
os.makedirs(output_dir, exist_ok=True)
stan_file = "./multipleLR.stan"
exe_file = "./multipleLR"
stanファイル
array型を使用した場合
data {
// 入力のチェックにもなるのでわかるところは下限と上限をつける
int<lower=0> N;
array[N] real<lower=0, upper=1> A;
array[N] real<lower=0, upper=1> Score;
array[N] real<lower=0, upper=1> Y;
}
parameters {
real b1;
real b2;
real b3;
real<lower=0> sigma;
}
transformed parameters {
array[N] real mu;
for (n in 1:N) {
mu[n] = b1 + b2*A[n] + b3*Score[n];
}
}
model {
// 各パラメータの無情報事前分布の指定は省略
for (n in 1:N) {
Y[n] ~ normal(mu[n], sigma);
}
}
generated quantities {
/* 各サンプルの予測分布を得るために、ここで正規分布から乱数を発生させて計算 */
array[N] real y_pred;
for (n in 1:N) {
y_pred[n] = normal_rng(mu[n], sigma);
}
}
matrix型を使用した場合
data {
// 入力のチェックにもなるのでわかるところは下限と上限をつける
int<lower=0> N;
int<lower=0> K;
matrix<lower=0, upper=1>[N, K] X;
vector<lower=0, upper=1>[N] y;
}
parameters {
real alpha;
vector[K] beta;
real<lower=0> sigma;
}
transformed parameters {
vector[N] mu;
mu = alpha + X*beta;
}
model {
// 各パラメータの無情報事前分布の指定は省略
y ~ normal(mu, sigma);
}
generated quantities {
/* 各サンプルの予測分布を得るために、ここで正規分布から乱数を発生させて計算 */
array[N] real y_pred;
y_pred = normal_rng(mu, sigma);
}
コンパイル
if not os.path.exists(exe_file):
model = CmdStanModel(stan_file=stan_file)
else:
model = CmdStanModel(exe_file=exe_file)
22:51:04 - cmdstanpy - INFO - compiling stan file /Users/yyamaguchi/Desktop/programming/BayesianStatistics/src/02_multipleLinearReg/multipleLR.stan to exe file /Users/yyamaguchi/Desktop/programming/BayesianStatistics/src/02_multipleLinearReg/multipleLR
22:51:08 - cmdstanpy - INFO - compiled model executable: /Users/yyamaguchi/Desktop/programming/BayesianStatistics/src/02_multipleLinearReg/multipleLR
データ
scaling
各パラメータのスケールがバラバラだとパラメータの探索が非効率になるため、Scoreのスケールを修正する
df["Score"] = df["Score"] / 200.
data = {
"N": len(df),
"A": df["A"].values,
"Score": df["Score"].values,
"Y": df["Y"].values,
}
print(data)
フィッティング
CmdStanPyでは以下のsample()メソッドにより、ハミルトンモンテカルロ(HMC)サンプリングを使って、データを条件としたモデルに対するベイズ推定を行います。このメソッドは、モデルとデータに対してStanのHMC-NUTSサンプラーを実行し、CmdStanMCMCオブジェクトを返します。
データは以下のようにPythonの辞書型でも渡せますし、ファイルパスを渡すこともできます。
import multiprocessing
num_cpu = multiprocessing.cpu_count()
fit = model.sample(
data=data,
chains=4, # chain数
seed=1, # seed固定
iter_warmup=1000, # warmupの数
iter_sampling=2000, # samplingの数
parallel_chains=num_cpu, # 並列数
save_warmup=True, # warmupもCSVに保存
thin=1, # サンプリング間隔
output_dir=output_dir, # 出力先
# show_console=True, # 標準出力
show_progress=True # progress出力
)
22:53:50 - cmdstanpy - INFO - Requested 8 parallel_chains but only 4 required, will run all chains in parallel.
22:53:51 - cmdstanpy - INFO - CmdStan start processing
chain 1 | | 00:00 Status
chain 2 | | 00:00 Status
chain 3 | | 00:00 Status
chain 4 | | 00:00 Status
22:53:51 - cmdstanpy - INFO - CmdStan done processing.
22:53:51 - cmdstanpy - WARNING - Non-fatal error during sampling:
Exception: normal_lpdf: Scale parameter is 0, but must be positive! (in '/Users/yyamaguchi/Desktop/programming/BayesianStatistics/src/02_multipleLinearReg/multipleLR.stan', line 26, column 8 to column 36)
Consider re-running with show_console=True if the above output is unclear!
結果の確認
fit.summary()で各パラメータの要約が見られます。lp__は対数事後確率で他のパラメータ同様に収束する必要があります。
また、N_EffはStanが自己相関等から判断した実効的なMCMCサンプル数です。この数が少ないと収束しづらいパラメータであることが分かります。大体1000サンプルほどあると良いものです。
fit.summary()
arvizによる収束性判断
import numpy as np
import arviz as az
import xarray as xr
az.style.use("arviz-darkgrid")
cmdstanpyからarviz用のデータへ変換
cmdstanpy_data = az.from_cmdstanpy(
posterior=fit,
log_likelihood="lp__",
)
cmdstanpy_data
ll_data = cmdstanpy_data.log_likelihood
対数事後確率のプロット
az.plot_trace(ll_data, compact=False);
自己相関プロット
ほとんど自己相関がないことがわかる
az.plot_autocorr(cmdstanpy_data, grid=(4, 4), var_names=["b1", "b2", "b3", "sigma"]);
事後分布のプロット
どのパラメータも95%信頼区間は0より大きい確率を持っているため、確かにb3が大きければYも大きいというように言うことができる。このようなパラメータの値が０より大きい確率をBayesian p-valueと呼ぶ。
az.plot_posterior(cmdstanpy_data, var_names=["b1", "b2", "b3", "sigma"]);
トレースプロット
ある点を中心にプロットされているので収束していることがわかる。
az.plot_trace(cmdstanpy_data, compact=False, var_names=["b1", "b2", "b3", "sigma"]);
パラメータのペアプロット
aとbは負の相関があることがわかる。
az.plot_pair(
cmdstanpy_data,
divergences=True,
var_names=["b1", "b2", "b3", "sigma"]
);
実測値と予測値プロット
df_summary = fit.summary()
df_summary = df_summary.loc[df_summary.index.str.contains("y_pred"), :]
df_summary["y_act"] = df["Y"].values
def plot_yy(y_act, y_pred, err_lower, err_upper):
min_value = min(y_act) - np.abs(min(y_act)*0.3)
max_value = max(y_act) + np.abs(max(y_act)*0.3)
plt.figure(figsize=(5, 5))
plt.plot([min_value, max_value], [min_value, max_value], c="lightgray")
plt.scatter(y_act, y_pred)
plt.errorbar(y_act, y_pred, yerr=[err_lower, err_upper], fmt="o")
plt.xlim([min_value, max_value])
plt.ylim([min_value, max_value])
plt.xlabel("y_actual")
plt.ylabel("y_pred")
plt.show()
plot_yy(y_act=df_summary["y_act"], y_pred=df_summary["Mean"], err_lower=df_summary["Mean"]-df_summary["5%"], err_upper=-df_summary["Mean"]+df_summary["95%"])
推定されたノイズの分布の確認
今回はAとScore以外の影響を正規分布に従うと想定した。そこで実際に
ここで、
結果としては大体同じような分布で想定していた分布に対し乖離は少なそうです。
df_summary = fit.summary()
df_summary = df_summary.loc[df_summary.index.str.contains("mu"), :]
array_epsilon = df["Y"].values - df_summary["50%"].values
from scipy.stats import norm
import statistics
plt.hist(array_epsilon, density=True)
x_axis = np.arange(-0.15, 0.15, 0.01)
plt.plot(x_axis, norm.pdf(x_axis, 0, fit.summary().loc["sigma", "Mean"]))
Discussion