PythonでBayes 編
Pythonを使ったベイズ統計
Yama-Lab.統計 ベイズ統計編を開始しました。
Google Coraboratory
Googleドライブから新規作成 → Google Coraboratory （"新規作成"の"その他"に"Google Coraboratory"がない場合は、"アプリを追加"で"Google Coraboratory"を追加しておく）
GoogleDriveのマウント方法 → Googleドライブとの連携について許可すると、フォルダアイコンの"content"に"drive"がでてくる。ファイルの保存や読み込みに利用する。
必要なライブラリのインポート
# 計算を素早くしてくれる
import numpy as np
# フレーム（エクセル）で行列の計算や集計など
import pandas as pd
# 機械学習
import seaborn as sns
# MCMC
import pystan
# グラフ
import matplotlib.pyplot as plt
# Pythonの中で図表示
%matplotlib inline
# 表示桁数の指定
%precision 3
# 複雑な統計計算
from scipy.stats import mstats
データの読み込み
data = pd.read.csv ("ファイル名.csv", sep=";", na_values=".") # csvファイルをdataという名前をつけて読み込む。第２引数以降はオプション
data.shape # データの数
pd.DataFrame.describe ( )
data.head ( ) # データの最初の５行目までを表示（きちんと読み込まれているか確認）
マックの場合で文字化けする場合
data = pd.read_csv("/content/drive/MyDrive/Python/フォルダ名/ファイル名.csv", encoding='utf-8')
相関
x = data ["age"] # ｘにdataの中の"age"の列を指定
y = data ["P1"] # y にdataの中の"P1"の列を指定
np.corrcoef(x, y) # Numpyで相関係数を確認
作図
散布図
plt.scatter(data ["P1"], ["age"]) #P1 と age の散布図
Stan 単回帰分析
# 単回帰分析 age でP2を回帰
stan_model = """
data {
int N;
real X [N];
real Y [N];
}
parameters {
real a;
real b;
real , low = 0 sigma;
}
model {
for (n in 1 1:N) {
Y [N] ~ normal (a * X[N] + b, sigma);
}
}
"""
sm = pystan.StanModel (model_code=stan_model)
stan_data = { "N" : df.shape[0], "X" : df["age"], "Y" : df["P2"] }
fit = sm.sampling(data = stan_data, iter=2000, warmup=500, chains = 3) # 繰り返し数2,000回、ウォームアップ500回、３セット繰り返し
fit
※ 確認する結果
１．Rhat値 1.0がベスト
各値を確認
fig = fit.plot() # 事後分布とMCMCの結果
結果の平均値から予測直線を引く
a = 0.16
b =-0.89
x = np.arange(65, 100, 1)
y = a * x + b
plt.plot(x, y)
plt.scatter(df["age"], df["P2"])
ベイズ信用区間
stan_model = """
data {
int N ;
real X [N];
real Y [N];
int N_s;
real X_s [N_s];
}
parameters {
real a;
real b;
real < lower = 0 > sigma;
}
model {
for (n in 1:N) {
Y [n] ~ normal (a * X[n] + b, sigma);
}
}
generated quantities {
real Y_s [N_s];
for (n in 1 : N_s) {
Y_s [n] = normal_rng (a * X_s [n] + b, sigma) ;
}
}
"""
sm = pystan.StanModel(model_code=stan_model)
X_s = np.arange( 65, 100, 1)
N_s = X_s.shape[0]
stan_data = { "N" : df.shape[0], "X" : df["P2"], "Y" : df["age"] , "N_s" : N_s, "X_s" : X_s}
fit = sm.sampling(data = stan_data, iter = 2000, warmup = 500, chains = 3)
fit.extract("a") # a についてMCMCで取得したデータを抽出
ms_a = fit.extract("a")["a"]
plt.hist(ms_a)
ms_b = fit.extract("b")["b"]
df_b = pd.DataFrame([])
for i in range(65, 100,1):
df_b[i] = ms_a * i + ms_b
df_b
low_y50, high_y50 = mstats.mquantiles(df_b, [0.25, 0.75], axis = 0)
low_y95, high_y95 = mstats.mquantiles(df_b, [0.025, 0.975], axis = 0)
plt.scatter(df["age"], df["P2"])
plt.fill_between(X_s, low_y50, high_y50, alpha = 0.6, color = "darkgray")
plt.fill_between(X_s, low_y95, high_y95, alpha = 0.3, color = "gray")
a = 0.16
b =-0.89
y = a * X_s + b
plt.plot(X_s, y, color = "black")
ベイズ予測区間
Y_p = fit.extract("Y_s")["Y_s"]
low_y, high_y = msats.mquantiles((Y_p, [0.025,0.975], asis = 0)
plt.scatter(df["age"], df["P2"])
plt.fill_between(X_s, low_y, high_y, alpha = 0.3, color = "gray")
a = 0.16
b =-0.89
y = a * X_s + b
plt.plot(X_s, y, color = "black")
重回帰分析
必要な変数のみのデータフレームを作成
dfm = df[["P1", "age", "economy", "sex"]]
変数間の関係を図で確認
g = sns.PairGrid (dfm)
g= g.map_lower (sns.kdeplot)
g = g.map_diag (sns.distplot, kde = False)
g = g.map_upper (plt.scatter)
stan_model = """
data {
int N ;
real age [N];
real economy [N];
real sex [N];
real Y [N];
}
parameters {
real a ;
real e;
real s;
real b;
real < lower = 0 > sigma;
}
model {
real mu;
for (n in 1:N) {
mu = a * age[n] + e * economy [n] + s * sex [n] + b;
Y [n] ~ normal (mu, sigma);
}
}
"""
sm = pystan.StanModel(model_code=stan_model)
stan_data = { "N" : df.shape[0], "age" : df ["age"], "economy" : df ["economy"], "sex" : df["sex"], "Y" : df["P1"] }
fit = sm.sampling(data = stan_data, iter = 2000, warmup = 500, chains = 4)
fit
fig = fit.plot()