Pythonでデータサイエンスを試す(8_MCMC)

前回から大分、間があいてしまいましたが、
Pythonでデータサイエンス続き編です。

前回のアソシエーションルールまでで、
本編は一通り終了しました。

残りは10章に載っている、
様々な手法に触れてみよう編です。

手を動かしながら学ぶ ビジネスに活かすデータマイニング

今回はベイジアンモデリング
教科書はこれ、ですが、ベイジアンモデリングそのものの説明は最小限でした。

準備

Pythonでベイジアンモデリングを用いるには、
MCMCを扱えるpystanを使用します。
これは重力波の研究にも使われたツールで、
StanというMCMCを扱うライブラリのPythonラッパーです。

andrewgelman.com

pipで入ります。

pip install pystan

そもそもベイジアンモデリングとMCMCとは

すごくざっくりと自分の理解を書いておくと、
GLMでは扱えない複雑なモデルを扱うための手法、です。
そのモデルの推定を行うために、行うシミュレーションがMCMCです。

ちゃんとした理論は、
専門書を読んできちんと学んだ方が良いと思います。

定番ですが、緑本。

データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (確率と情報の科学)

下準備

import pysta
import pandas as pd
import matplotlib.pyplot as plt

今回データは緑本の著者さんのHPよりとってくる、
とのことなので、wgetでとってきます。
（普通にダウンロードしてもいいです）

wget http://hosho.ees.hokudai.ac.jp/~kubo/stat/iwanamibook/fig/hbm/data7a.csv

dat = pd.read_csv("data7a.csv")
dat.head()

 id  y
0   1   0
1   2   2
2   3   7

使ってみる

pystanの使い方は、
stanのインターフェイスなので、
stan用のファイルと、
そこにデータを流し込んで、キックするpythonスクリプトを作ります。

つまり、stanわかんと扱えないということか・・・

stanファイルは本書のRstan用と同じであることを祈って、
本書に書かれてるとおりに作成しました。

data{
    int<lower=0>N;
    int<lower=0>y[N];
}

parameters{
    real beta;
    real r[N];
    real<lower=0>s;
}

transformed parameters{
    real q[N];
    for(i in 1:N)
        q[i] <- inv_logit(beta+r[i]);
}

model{
    for (i in 1:N)
        y[i] ~ binomial(8,q[i]);
    beta ~ normal(0,1.0+2);
    for(i in 1:N)
        r[i]~normal(0,s);
    s~uniform(0,1.0e+4);

}

シミュレーションを行うための
基データと変動させるパラメータを定義。
その後パラメータの変動のさせ方を定義しているっぽい。
ベータ分布と、ロジット関数を使って、最後のmodelでモデルの最適化をしてるんでしょうか。
（コード読んだ妄想なので、多分間違ってる)

で、データをインターフェイスするコードを書いてみます。
pystanの説明を参考に書いてみます。

N = len(df["y"]) #StanファイルのN用のデータ
y = df["y"].values #Stanファイルのy用のデータ

stan_data = {"N":N,"y":y} #pythonのデータとstanファイルの変数を紐付ける

# そしてキック
fit = pystan.stan(
    file='Stanファイル', # 使用するstanファイル
    data=stan_data,# 上で定義したデータ対応をdataとして渡す
    iter=1000,# 試行回数
    chains=4 # 連鎖回数
)

fit.plot()
fit.summary()

OrderedDict([('c_summary',
              array([[[  5.91655341e-02,   3.33755865e-01,  -6.68194488e-02,
                         3.54550700e-01],
                      [  6.44419330e-02,   3.60163156e-01,   7.64982789e-02,
                         3.18692525e-01],
                      [ -6.23284252e-01,  -1.67281761e-01,   5.92567233e-02,
                         2.87728466e-01],
   mcmc                   ...,

動いた!
けど、読み取り方がいまいちわからないので、
MCMCとstanの勉強をします。。。
（ちなみに結構マシンパワーを食いました)

残りはanova、時系列分析です。
10章は一つ一つが専門的だったりするので難しい。

前回はこちら

esu-ko.hatenablog.com