PyMC3で変動係数一定な回帰ベイズモデルを作ってみた

線形回帰の誤差分布

線形回帰分析を行うときには、一般的に誤差は平均ゼロ、標準偏差一定の正規分布に従う事を仮定します。
この仮定の下で、最尤法と最小二乗法が一致する事が知られているからです。

一方で、そのような仮定が明らかに成り立たない場合もあります。
例えば年齢で収入を予測するモデルを作るとしましょう。
この時、若年代のばらつきは小さく、老年代のばらつきは大きい事が明らかです。
なぜなら若くして高収入な人は殆ど居らず、一方で年寄りには金持ちも金持ちもいれば貧困にあえぐ層も多くいるからです。
このような場合、年代に依らず誤差分布が一定だと仮定するのは無理があることがわかります。
一方で、変動係数は年代によらず一定だと仮定しても良いかもしれません。
このようなモデルを考えてみます。

変動係数

変動係数cvとは、標準偏差を平均値で割った統計量のことです。
$ cv =\frac{\sigma}{\bar{x}} $
この統計量は無次元なので、単位によらずばらつきの比較をするときに使えます。

一般化線形モデルとベイズ統計モデル

このように、誤差が(一定の)正規分布に従わない場合、伝統的な統計解析では、一般化線形モデルという枠組みを考えることが一般的です。
しかしこの一般化線形モデルには制限があり、今回のベイズ解析のように直接的に任意の誤差分布を指定する事は難しいものになっています。
又、一般化線形モデルでは「誤差の大きさそのもの」を推定することは困難です。(理論的にはやろうと思えば出来そうですが、そこまでサポートしているライブラリが見つかりませんでした)
したがって、今回はより一般的な手法である、ベイズ統計をベースにしたPyMC3を使用してモデルを構築してみました。

変動係数一定のデータを作成

xが大きくなるにしたがって、ばらつきが大きくなっていくデータを作ります
jupyterで適当の次のコードを実行してください。

data.ipynb

import pymc3 as pm
import numpy as np
import pandas as pd
import scipy.stats as stats
import matplotlib.pyplot as plt

np.random.seed(314)
N=1000
beta_real=0.9
eps_real = np.random.normal(0, 0.3, size=N)

x=np.random.normal(3, 1, N)
y_real = beta_real * x
y = y_real + eps_real * x

plt.plot(x,y,"b.")
plt.xlabel("$x$", fontsize=16)
plt.ylabel("$y$", fontsize=16, rotation=0)
plt.plot(x, y_real, "k")

実行結果

このように、xの大きさに従ってばらつきが大きくなっていくデータが作成できました。

#変動係数一定なベイズモデルを生成
上で作ったデータを予測するベイズモデルを作りました。

model.ipynb

with pm.Model()as model:
    beta = pm.Normal("beta", mu=0, sd=1) # 回帰係数の事前分布
    cv = pm.HalfCauchy("cv", 1) # 変動係数の事前分布
    mu = pm.Deterministic("mu", beta*x) # 回帰モデル。線形なモデルを仮定する。
    # 誤差分布モデル。変動係数cv = sd/|x|が一定なモデルを生成する。
    y_pred = pm.Normal("y_pred", mu=mu, sd=np.abs(x)*cv, observed=y)
    
    #以下ボイラーコード
    start = pm.find_MAP()
    step = pm.Metropolis()
    trace = pm.sample(11000, step, start, njobs=1)

trace_n = trace[100:]
pm.traceplot(trace_n)

実行結果

グラフはパラメータの事後分布です。真の値であるbeta=0.9, cv=0.3に大体ピークがきてますね。モデル構築は上手くいっているみたいです。

予測区間

単なる線形回帰モデルを作るだけではベイズを使ったうまみがありません。
次に区間予測をしてみましょう

interval_reg.ipynb

idx = np.argsort(x)
x_ord=x[idx]
beta_m=trace_n["beta"].mean()
plt.plot(x,y,"b.")
plt.plot(x, beta_m * x, c="k", label="y={:.2f}*x".format(beta_m))
sig0 = pm.hpd(ppc["y_pred"], alpha=0.5)[idx] # 50パーセント信頼区間予測
sig1 = pm.hpd(ppc["y_pred"], alpha=0.05)[idx] # 95パーセント信頼区間予測
plt.fill_between(x_ord, sig0[:,0], sig0[:,1], color="gray", alpha=1) 
plt.fill_between(x_ord, sig1[:,0], sig1[:,1], color="gray", alpha=0.5)

plt.xlabel("$x$", fontsize=16)
plt.ylabel("$y$", fontsize=16, rotation=0)
plt.legend()

実行結果

グラフの濃いグレーは50%, 薄いグレーは95%信頼区間を表しています。
xに比例して区間幅が広くなっていることがわかりますね。
データにもうまいこと当てはまっていることが確認できます。

参考文献

Pythonによるベイズ統計モデリング: PyMCでのデータ分析実践ガイド - Osvaldo Martin (原著), オズワルド マーティン (著), 金子 武久 (翻訳)