「Pythonによるベイズ統計モデリング」まとめ #3

「Pythonによるベイズ統計モデリング―PyMCでのデータ分析実践ガイド― 」：

Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド

前回↓ 

ryosuke-okubo.hatenablog.com

#3では第4章と第5章を扱う。データは一貫してseabornのirisを用いる。

下準備

import pymc3 as pm
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('seaborn-darkgrid')
np.set_printoptions(precision=2)
pd.set_option('display.precision', 2)

 データ

iris = sns.load_dataset('iris')
iris.head()

 

sns.pairplot(iris, hue='species', diag_kind='kde') plt.figure()

第4章　線形回帰モデルによるデータの理解と予測

線形単回帰

petal_widthとpetal_lengthの値から線形単回帰を行ってみる。

df = iris.query("species == ('setosa', 'versicolor')")
x_n = 'petal_width' 
x = df[x_n].values
y_n = 'petal_length'
y = df[y_n].values
plt.scatter(x, y)
plt.xlabel('$width$', fontsize=16)
plt.ylabel('$length$', fontsize=16)

plt.figure()

線形モデルの式は

で表される。これは，確率的には次のように表現される。

今回の例では，各パラメータを次のように設定する。

，

 

with pm.Model() as model:
    alpha = pm.Normal('alpha', mu=0, sd=1)
    beta = pm.Normal('beta', mu=0, sd=1)
    epsilon = pm.HalfCauchy('epsilon', 5)

    mu = pm.Deterministic('mu', alpha + beta * x)
    y_pred = pm.Normal('y_pred', mu=mu, sd=epsilon, observed=y)

    start = pm.find_MAP()
    step = pm.Metropolis()
    trace = pm.sample(11000, step, start, njobs=1)
    
trace_n = trace[1000:]
pm.traceplot(trace_n)

plt.figure()

線形回帰では，しばしば激しい自己相関がみられる。解決法として，

とデータベクトルを中心化することがあげられる。この際，各パラメータは

，

と変換されている。

事後分布の要約

ppc = pm.sample_ppc(trace_n, samples=1000, model=model) #事後予測サンプルを得る

alpha_m = trace_n['alpha'].mean()
beta_m = trace_n['beta'].mean()
plt.plot(x, y, 'b.');
idx = np.argsort(x)
plt.plot(x, alpha_m + beta_m * x, c='k', label='y = {:.2f} + {:.2f} * x'.format(alpha_m, beta_m))

sig0 = pm.hpd(ppc['y_pred'], alpha=0.5)[idx]
sig1 = pm.hpd(ppc['y_pred'], alpha=0.05)[idx]
plt.fill_between(x_ord, sig0[:,0], sig0[:,1], color='gray', alpha=1) #50%HPD
plt.fill_between(x_ord, sig1[:,0], sig1[:,1], color='gray', alpha=0.5) #95%HPD

plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16, rotation=0)
plt.savefig('img410.png', dpi=300, figsize=(5.5, 5.5))

plt.figure()

その他，本書には次の項目がある。

• 頑健線形回帰
• 階層線形回帰
• 多項式回帰
• 線形重回帰 

第5章　ロジスティック回帰による結果変数の分類

 ロジスティック回帰

ロジスティックモデルは次のように記述される。

ここで

はシグモイド関数である。

z = np.linspace(-10, 10, 100)
logistic = 1 / (1 + np.exp(-z))
plt.plot(z, logistic)
plt.xlabel('$z$', fontsize=18)
plt.ylabel('$logistic(z)$', fontsize=18)

plt.figure()

例として，petal_lengthの値からsetosaとversicolorを分けることにする。

df = iris.query("species == ('setosa', 'versicolor')")
y = pd.Categorical(df['species']).codes
x_n = 'petal_length' 
x = df[x_n].values
plt.xlabel(x_n, fontsize=16)
plt.ylabel(r'$\theta$', rotation=0, fontsize=16)

plt.scatter(x, y)

モデル作成

with pm.Model() as model_0:
    alpha = pm.Normal('alpha', mu=0, sd=1)
    beta = pm.Normal('beta', mu=0, sd=1)
  
    mu = alpha + pm.math.dot(x, beta)
    theta = pm.Deterministic('theta', 1 / (1 + pm.math.exp(-mu)))
    bd = pm.Deterministic('bd', -alpha/beta)
    yl = pm.Bernoulli('yl', p=theta, observed=y)

    trace = pm.sample(5000)

chain = trace[1000:]
varnames = ['alpha', 'beta', 'bd']
pm.traceplot(chain, varnames)

plt.figure()

pm.summary(chain, varnames)

推定結果

theta = chain['theta'].mean(axis=0)
idx = np.argsort(x)
plt.plot(x[idx], theta[idx], color='b', lw=3);
plt.axvline(chain['bd'].mean(), ymax=1, color='r')
bd_hpd = pm.hpd(chain_0['bd'])
plt.fill_betweenx([0, 1], bd_hpd[0], bd_hpd[1], color='r', alpha=0.5)

plt.plot(x, y, 'o', color='k')
theta_hpd = pm.hpd(chain['theta'])[idx]
plt.fill_between(x[idx], theta_hpd[:,0], theta_hpd[:,1], color='b', alpha=0.5)

plt.xlabel(x_n, fontsize=16)
plt.ylabel(r'$\theta$', rotation=0, fontsize=16)

plt.figure() 

 その他，本書には次の項目がある。

•  多重ロジスティック回帰
• ソフトマックス回帰
• 線形判別分析

次回↓

ryosuke-okubo.hatenablog.com