PyStanでOLS①
今回は、流行りのStanをPythonから使ってみようと思う。PyStanを使うことでPythonからStanが呼べるようになる。
Stanは、ハミルトニアンモンテカルロ(HMC)とメトロポリス・ヘイスティングス(MH)の合わせ技で、従来のMHアルゴリズムによるMCMCより効率的なサンプリングが可能と言われている。
さらに、NUTS(No-U-Turn Sampler)という技が使われていて、これはHMCを実装する際にユーザーが決めるべきパラメータを自動で決めてくれる優れもの。
この辺は、理論的裏付けが下記の論文に書いてあるようなので、いつか真面目に勉強したい。
http://www.stat.columbia.edu/~gelman/research/published/nuts.pdf
Stanは、RやPython、Matlabから呼べて、RからStanを呼ぶRStanはそれなりに触ったことがあるが、PyStanは使ったことがないので、今回触ってみる。
推定に使うモデルはいつもの超簡単な単回帰式。
推定したいモデルは、Pythonコードとは別に用意したStanコードの方に書いておく。
Pythonコードからは、Stanコードにデータを渡してキックするイメージ。
PyStanを実行することで得られるパラメータ推定値の分布とパスのプロットがこちら。
上段が定数項と回帰係数、下段が残差の標準偏差となっている。定数項と回帰係数はベクトルにまとめたため二つ同時に出ている。
続いて、パラメータ推定値のサマリー。
95%信用区間に真のパラメータ（ベイズの考え方からは本来こういう言い方はしないが‥）が含まれていることがわかる。
In [105]: fitchan
Out[105]:
Inference for Stan model: anon_model_b665a1feb6ebc5ee8386c19e8517f522.
1 chains, each with iter=2000; warmup=100; thin=1;
post-warmup draws per chain=1900, total post-warmup draws=1900.
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
beta[0] 2.04 2.3e-3 0.04 1.96 2.01 2.04 2.07 2.13 352.0 1.0
beta[1] 1.47 1.6e-3 0.03 1.41 1.45 1.47 1.49 1.53 381.0 1.0
sigma 0.97 9.5e-4 0.02 0.93 0.95 0.97 0.99 1.01 555.0 1.0
lp__ -468.7 0.07 1.25 -471.8 -469.2 -468.3 -467.8 -467.3 312.0 1.0
Samples were drawn using NUTS(diag_e) at 02/27/15 17:19:26.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).一応、OLSの結果とも比較しておくと、ほぼ完全に一致していることが確認できる。
Stanのパラメータの平均値と標準偏差は、OLSのパラメータ推定値と標準誤差にそれぞれ相当する。
print(results.summary())
OLS Regression Results
==============================================================================
Dep. Variable: y R-squared: 0.692
Model: OLS Adj. R-squared: 0.692
Method: Least Squares F-statistic: 2241.
Date: Fri, 27 Feb 2015 Prob (F-statistic): 2.28e-257
Time: 17:27:16 Log-Likelihood: -1386.1
No. Observations: 1000 AIC: 2776.
Df Residuals: 998 BIC: 2786.
Df Model: 1
==============================================================================
coef std err t P>|t| [95.0% Conf. Int.]
------------------------------------------------------------------------------
const 2.0431 0.043 47.942 0.000 1.959 2.127
x1 1.4691 0.031 47.339 0.000 1.408 1.530
==============================================================================
Omnibus: 1.076 Durbin-Watson: 2.074
Prob(Omnibus): 0.584 Jarque-Bera (JB): 1.155
Skew: 0.059 Prob(JB): 0.561
Kurtosis: 2.883 Cond. No. 2.53
==============================================================================Pythonコードはこちら。
# -*- coding: utf-8 -*- from __future__ import print_function, division import matplotlib.pyplot as plt import numpy as np import seaborn as sns import statsmodels.api as sm import pystan #サンプリング数 N=1000 #データ生成 np.random.seed(0) x=np.random.normal(1,1,N) eps=np.random.normal(0,1,N) y=2+1.5*x+eps #OLS(statsmodels) X=sm.add_constant(x, prepend=True) model = sm.OLS(y,X) results = model.fit() print(results.summary()) #Stan N = X.shape[0] K = X.shape[1] stan_data = {'N': N, 'K': K, 'y': y, 'X': X} model = pystan.StanModel('D:\Python/OLS.stan') fitchan = model.sampling(data=stan_data, iter=2000, warmup=100,chains=1) fitchan.plot() fitchan
Stanコードはこちら。
data {
int<lower=1> N;
int<lower=1> K;
real y[N];
matrix[N,K] X;
}
parameters {
vector[K] beta;
real<lower=0.0001> sigma;
}
model {
y ~ normal(X*beta,sigma);
}
