前回 : Pythonで学ぶ！ベイズ統計第7回 - ma-hi-ro Hobby

Pythonで学ぶ！ベイズ統計第8回

前回までは、確率分布のパラメータを推定していた
今回は、データに対して、線形モデルのパラメータ推定を MCMC で行う

Pythonで学ぶ！ベイズ統計第8回

データの生成

自己相関の検証のため、平均のズレたデータ、で、
ロバスト線形回帰の検証のため、外れ値のあるデータ、にする

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(5, 15, 30)

np.random.seed(1)
err = np.random.normal(0, 1, size=30)

y_true = 1 + 2 * x
y = y_true + err
y[-2] = 50
y[-4] = 60

fig, ax = plt.subplots()
ax.scatter(x, y, label='data')
ax.plot(x, y_true, 'k--', label='true')
ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.legend()
ax.grid()
fig.tight_layout()

$y = 1 + 2x$ に $\sigma = 1$ のノイズが追加され、外れ値を２つ含んでいる
$x$ の平均は 10 になっていて、平均値を除去しないと、自己相関が見えてくるはず

線形回帰（そのまま）

MCMCで回帰を行う
データの分布は正規分布を設定する

import pymc3 as pm
import warnings
warnings.simplefilter('ignore', FutureWarning)

with pm.Model() as model:
    beta0 = pm.Normal('beta0', mu=0, sigma=100)
    beta1 = pm.Normal('beta1', mu=0, sigma=100)
    epsilon = pm.HalfNormal('epsilon', sigma=5)
    mu = pm.Deterministic('mu', beta0 + beta1 * x)

    pm.Normal('y', mu=mu, sigma=epsilon, observed=y)
    trace = pm.sample(random_seed=1)

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 243 seconds.

names = ['beta0', 'beta1', 'epsilon']
pm.plot_trace(trace, var_names=names)

pm.summary(trace, var_names=names)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta0	-5.500	4.377	-13.779	2.723	0.121	0.094	1335.0	1277.0	1.00
beta1	2.815	0.420	2.001	3.572	0.012	0.008	1330.0	1272.0	1.00
epsilon	6.384	0.883	4.936	8.185	0.029	0.021	815.0	438.0	1.01

切片（beta0）の精度が悪い
自己相関を見てみると、下記のように相関が大きく、MCMCができていないことがわかる

pm.plot_autocorr(trace, var_names=names, combined=True)

線形回帰（平均値除去）

$x$ の平均値を除去する

x_c = x - x.mean()
with pm.Model() as model_c:
    beta0 = pm.Normal('beta0', mu=0, sigma=100)
    beta1 = pm.Normal('beta1', mu=0, sigma=100)
    epsilon = pm.HalfNormal('epsilon', sigma=5)
    mu = pm.Deterministic('mu', beta0 + beta1 * x_c)

    pm.Normal('y', mu=mu, sigma=epsilon, observed=y)
    trace_c = pm.sample(random_seed=1)

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 143 seconds.

names = ['beta0', 'beta1', 'epsilon']
pm.plot_trace(trace_c, var_names=names)

pm.summary(trace_c, var_names=names)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta0	22.629	1.163	20.481	24.805	0.018	0.012	4388.0	2904.0	1.0
beta1	2.830	0.392	2.117	3.585	0.006	0.004	4805.0	2988.0	1.0
epsilon	6.387	0.859	4.905	7.979	0.014 0.010	3919.0	2534.0	1.0

今度は、切片（beta0）の精度が改善している

自己相関は下図

pm.plot_autocorr(trace_c, var_names=names, combined=True)

本来の切片の値は下記で求まる

beta0 = trace_c['beta0'] - trace_c['beta1'] * x.mean()
beta0.mean()

-5.671645097099628

MCMCでは、信頼区間を求めることができる

with model_c:
    pp = pm.sample_posterior_predictive(trace_c)

fig, ax = plt.subplots()

ax.scatter(x, y, label='data')
ax.plot(x, trace_c['mu'].mean(axis=0), 'k', label='mean')

pm.plot_hdi(x, pp['y'], ax=ax, fill_kwargs={'color': 'gray', 'label': r'$94\%$ HDI', 'alpha': 0.2})
pm.plot_hdi(x, pp['y'], ax=ax, fill_kwargs={'color': 'gray', 'label': r'$50\%$ HDI', 'alpha': 0.5}, hdi_prob=0.5)

ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.legend()
ax.grid()
fig.tight_layout()

外れ値に引っ張られて精度が悪くなっている

ロバスト線形回帰

今度は、データの分布を「スチューデントのt分布」で与える
また、t分布のパラメータnuは指数分布で与える

x_c = x - x.mean()
with pm.Model() as model_t:
    beta0 = pm.Normal('beta0', mu=0, sigma=100)
    beta1 = pm.Normal('beta1', mu=0, sigma=100)
    nu = pm.Exponential('nu', 1/30)
    sigma = pm.HalfNormal('sigma', 5)
    mu = pm.Deterministic('mu', beta0 + beta1 * x_c)

    pm.StudentT('y', mu=mu, sigma=sigma, nu=nu, observed=y)
    trace_t = pm.sample(random_seed=1)

Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 281 seconds.

pm.summary(trace_t, var_names='~mu')

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
beta0	20.953	0.233	20.551	21.425	0.004	0.003	4204.0	2508.0	1.0
beta1	2.078	0.087	1.912	2.236	0.002	0.001	2505.0	2486.0	1.0
nu	1.437	0.497	0.625	2.395	0.010	0.007	2211.0	1980.0	1.0
sigma	0.916	0.252	0.457	1.390	0.005	0.004	1949.0	1471.0	1.0

beta0 = trace_t['beta0'] - trace_t['beta1'] * x.mean()
beta0.mean()

0.1772918671759852

beta0 と beta1 の推定も推定できている
（beat0が多少ズレているが、外れ値があるため、このぐらいか・・・）

with model_t:
    pp_t = pm.sample_posterior_predictive(trace_t)

fig, ax = plt.subplots()

ax.scatter(x, y, label='data')
ax.plot(x, trace_t['mu'].mean(axis=0), 'k', label='mean')

pm.plot_hdi(x, pp_t['y'], ax=ax, fill_kwargs={'color': 'gray', 'label': r'$94\%$ HDI', 'alpha': 0.2})
pm.plot_hdi(x, pp_t['y'], ax=ax, fill_kwargs={'color': 'gray', 'label': r'$50\%$ HDI', 'alpha': 0.5}, hdi_prob=0.5)

ax.set_xlabel(r'$x$')
ax.set_ylabel(r'$y$')
ax.legend()
ax.grid()
fig.tight_layout()