書記が数学やるだけ#782 多変量正規分布モデル・線形回帰モデルのベイズ統計
2変数データのベイズ統計について，多変量正規分布モデルと線形回帰モデルを実装してみる。
問題
回帰分析では恒例の樹木データを用いる：
説明
多変量正規分布について：
多次元正規モデルにおいて，事後分布は逆ウィシャート分布である。
ウィシャート分布は複雑な形をしており，これはカイ二乗分布の拡張である。
次に線形回帰モデル，これは事後分布は正規分布で表される。
多次元正規分布モデルと線形回帰モデルは，条件付き期待値の観点から密接に結びついている：
解答
import pymc3 as pm
import arviz as az
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')
#データの読み込み，予め用意しているものとする
data = pd.read_csv('回帰分析-1.csv')
data一旦は最尤推定をしておく。
import statsmodels
import statsmodels.api as sm
import statsmodels.formula.api as smf
#最尤推定
model = smf.ols('VOLUME ~ HEIGHT', data=data).fit()
model.summary()
最尤推定による線形回帰モデルを図示しておく。
x = data["HEIGHT"]
y = data["VOLUME"]
y_ML = -87.1236 + 1.5433 * x
plt.figure(figsize=(10, 5))
plt.subplot(1, 2, 1)
plt.plot(x, y, 'b.')
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16, rotation=0)
plt.plot(x, y_ML, 'k')
plt.subplot(1, 2, 2)
sns.kdeplot(y)
plt.xlabel('$y$', fontsize=16)
plt.figure()
ベイズ推定の前処理として標準化（平均0，分散1）を行なっておく。
data_s = stats.zscore(data)
data_s先に線形回帰モデルから構築する。モデルとして，事前分布を経験的に定める構造を採用する。
with pm.Model() as model_t:
alpha = pm.Normal('alpha', mu=0, sd=100)
beta = pm.Normal('beta', mu=0, sd=1)
epsilon = pm.HalfCauchy('epsilon', 5)
nu = pm.Deterministic('nu', pm.Exponential('nu_', 1/29) + 1)
y_pred = pm.StudentT('y_pred', mu=alpha + beta * x, sd=epsilon, nu=nu, observed=y)
start = pm.find_MAP()
step = pm.NUTS(scaling=start)
trace_t = pm.sample(1000, step=step, start=start)
az.plot_trace(trace)
plt.figure()
事後分布は以下の通り。
plt.plot(x, y, 'b.');
idx = range(0, len(trace['alpha']), 10)
plt.plot(x, trace['alpha'][idx] + trace['beta'][idx] * x[:,np.newaxis], c='gray', alpha=0.5);
plt.plot(x, alpha_m + beta_m * x, c='k', label='y = {:.2f} + {:.2f} * x'.format(alpha_m, beta_m))
plt.xlabel('$x$', fontsize=16)
plt.ylabel('$y$', fontsize=16, rotation=0)
plt.legend(loc=2, fontsize=14)
plt.figure()
続いて2変量正規モデルについて。
まず全体の観察として，以下の記事のコードを借りて表示した：
ではモデルを組んでいく。
with pm.Model() as pearson_model:
mu = pm.Normal('mu', mu=0, sd=1, shape=2)
sigma_1 = pm.HalfNormal('simga_1', 1)
sigma_2 = pm.HalfNormal('sigma_2', 1)
rho = pm.Uniform('rho', -1, 1)
cov = pm.math.stack(([sigma_12, sigma_1sigma_2rho], [sigma_1sigma_2rho, sigma_22]))
y_pred = pm.MvNormal('y_pred', mu=mu, cov=cov, observed=data)
start = pm.find_MAP()
step = pm.NUTS(scaling=start)
trace_p = pm.sample(10000, step=step, start=start)
pm.traceplot(trace_p)
plt.figure()
本記事のもくじはこちら：
