Pythonだけでベイジアンモデリングをしたい(pymc_線形回帰)

今回は線形回帰をモデリングしてみます。

係数と誤差はそれぞれ正規分布を過程します
これらを線形に計算する場合はDeterministicで確定的であることを宣言します。

これがさらに正規分布し、sdはコーシー分布を設定しておきます。

やりかた

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

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

やってみる

データの生成

import numpy as np
import pymc3 as pm

N=100
alpha_real = 2.5
beta_real = 0.9
error_real = np.random.normal(0,0.5,size=N)

x = np.random.normal(10,1,size=N)
y_real = alpha_real + beta_real * x

y = y_real + error_real

生成したデータを可視化してみます。

import matplotlib.pyplot as plt
plt.scatter(x,y)

モデリングと結果の確認

今回はサンプリングを自己相関に強いNUTSを用います。

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

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

with reg_mdl:
  start=pm.find_MAP()
  step=pm.NUTS()
  trace=pm.sample(1000,step=step,start=start,cores=4)

burn_n = 300
trace_n = trace[burn_n:]
pm.traceplot(trace_n)

最終結果の確認

pm.summary(trace_n,var_names=['alpha','beta','rand'])

 mean    sd  hdi_3%  hdi_97% mcse_mean   mcse_sd ess_mean    ess_sd  ess_bulk    ess_tail    r_hat
alpha   2.247   0.529   1.249   3.201   0.019   0.013   801.0   801.0   815.0   816.0   1.0
beta    0.919   0.052   0.823   1.015   0.002   0.001   797.0   785.0   811.0   795.0   1.0
rand    0.518   0.037   0.451   0.585   0.001   0.001   1400.0  1384.0  1426.0  1343.0  1.0

それっぽい結果がでました。r_hatも問題なさそうです。