Tech Tips: 3月 2017

2017年3月31日金曜日

PyStan入門（1）とりあえず動かす

　緑本を読んで、MCMCで統計モデルのパラメータ推定をすることに興味がわいてきたので、PyStanをさわっていこうと思う。
　まずは、とりあえずインストールして動かしてみました。

PyStanのインストール
pipで入ります。

$ pip install pystan

モデルの記述
stanのモデルは.stanという拡張子のファイルに書くのが一般的なようです。他にもPythonのコードにベタ書きするという方法もあるようです。
ここではstanファイルにモデルを記述することにし、以下の内容を8schools.stanというファイルに記述します。

data {
  int<lower=0> J; // number of schools 
  real y[J]; // estimated treatment effects
  real<lower=0> sigma[J]; // s.e. of effect estimates 
}
parameters {
  real mu; 
  real<lower=0> tau;
  real eta[J];
}
transformed parameters {
  real theta[J];
  for (j in 1:J)
    theta[j] <- mu + tau * eta[j];
}
model {
  eta ~ normal(0, 1);
  y ~ normal(theta, sigma);
}

上のモデルを以下のPythonコードから参照します。以下のコードをsample.pyという名前のファイルに記述します。

import pystan
import matplotlib.pyplot as plt

schools_dat = {
 'J': 8,
 'y': [28,  8, -3,  7, -1,  1, 18, 12],
 'sigma': [15, 10, 16, 11,  9, 11, 10, 18]
}

fit = pystan.stan(file='8schools.stan', data=schools_dat, iter=1000, chains=4)

print(fit)
fit.plot()
plt.show()

実行結果
実行してみます。

$ python sample.py
DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal): assignment operator <- deprecated in the Stan language; use = instead.

INFO:pystan:COMPILING THE C++ CODE FOR MODEL anon_model_286b3180dfa752c4cfedaf0241add0e4 NOW.
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 1)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 3)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 2)
Iteration:   1 / 1000 [  0%]  (Warmup) (Chain 0)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 1)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 2)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 3)
Iteration: 100 / 1000 [ 10%]  (Warmup) (Chain 0)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 1)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 0)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 0)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 1)
Iteration: 200 / 1000 [ 20%]  (Warmup) (Chain 2)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 1)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 0)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 3)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 3)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 3)
Iteration: 300 / 1000 [ 30%]  (Warmup) (Chain 2)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 0)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 1)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 0)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 1)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 3)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 0)
Iteration: 400 / 1000 [ 40%]  (Warmup) (Chain 2)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 1)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 3)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 1)
Iteration: 500 / 1000 [ 50%]  (Warmup) (Chain 2)
Iteration: 501 / 1000 [ 50%]  (Sampling) (Chain 2)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 0)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 3)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 1)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 3)
Iteration: 600 / 1000 [ 60%]  (Sampling) (Chain 2)
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 0)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 3)
# 
#  Elapsed Time: 0.044725 seconds (Warm-up)
#                0.037462 seconds (Sampling)
#                0.082187 seconds (Total)
# 
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 1)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 0)
Iteration: 700 / 1000 [ 70%]  (Sampling) (Chain 2)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 1)
# 
#  Elapsed Time: 0.05144 seconds (Warm-up)
#                0.040057 seconds (Sampling)
#                0.091497 seconds (Total)
# 
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 0)
# 
#  Elapsed Time: 0.04517 seconds (Warm-up)
#                0.046124 seconds (Sampling)
#                0.091294 seconds (Total)
# 
Iteration: 800 / 1000 [ 80%]  (Sampling) (Chain 2)
Iteration: 900 / 1000 [ 90%]  (Sampling) (Chain 2)
Iteration: 1000 / 1000 [100%]  (Sampling) (Chain 2)
# 
#  Elapsed Time: 0.053818 seconds (Warm-up)
#                0.040396 seconds (Sampling)
#                0.094214 seconds (Total)
# 
Inference for Stan model: anon_model_286b3180dfa752c4cfedaf0241add0e4.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

           mean se_mean     sd   2.5%    25%    50%    75%  97.5%  n_eff   Rhat
mu         7.81    0.17   5.09  -2.41   4.61   7.77  11.13  17.63    907    1.0
tau        6.63    0.24    5.6   0.29   2.45   5.33   9.17  21.19    531    1.0
eta[0]     0.42    0.02   0.94  -1.53  -0.17   0.44   1.05   2.17   1701    1.0
eta[1]     0.01    0.02   0.87  -1.67  -0.55 3.2e-3    0.6   1.74   1751    1.0
eta[2]    -0.23    0.02   0.89  -1.95  -0.81  -0.22   0.38   1.51   1866    1.0
eta[3]    -0.01    0.02    0.9  -1.82  -0.61-8.1e-3   0.58   1.79   1669    1.0
eta[4]    -0.33    0.02   0.87   -2.1  -0.87  -0.34   0.21   1.44   1603    1.0
eta[5]    -0.21    0.02    0.9   -1.9  -0.82  -0.23   0.39   1.55   1599    1.0
eta[6]     0.34    0.02    0.9  -1.42  -0.24   0.36   0.95   2.07   1537    1.0
eta[7]     0.06    0.02   0.94  -1.84  -0.56   0.07   0.66   1.87   1604    1.0
theta[0]  11.61    0.24   8.19  -1.43   6.15  10.55  15.65  31.15   1211    1.0
theta[1]   7.84    0.14   6.35  -5.56   4.05   7.76   11.7  20.83   1988    1.0
theta[2]   5.97     0.2    7.7  -11.7   1.85   6.38  10.91  20.08   1505    1.0
theta[3]   7.61    0.15   6.45  -5.51   3.46   7.67  11.86  20.35   1897    1.0
theta[4]   5.13    0.14   6.45  -8.65   0.99   5.58   9.56   17.0   2000    1.0
theta[5]   6.21    0.17   6.51  -7.34   2.29   6.46  10.55  18.69   1522    1.0
theta[6]  10.64    0.17   6.88  -1.89   6.09  10.06  14.68  26.06   1623    1.0
theta[7]   8.33     0.2   7.94  -7.88   3.77   8.07  12.54  25.84   1522    1.0
lp__      -4.82    0.11   2.66 -10.43  -6.54  -4.59  -2.84  -0.42    598    1.0

Samples were drawn using NUTS at Sat Apr  1 00:55:24 2017.
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).

C++のコンパイラが走って、MCMCのサンプリングが走って、パラメータの分布が出たようです。続いて各パラメータの分布図とiterationごとの変化がプロットされます。

細かいことは置いといて、何やら楽しそうなことが出来るということは分かりました。
次回は自分で簡単な乱数生成モデルを作って、それをモデリングしてパラメータ推定ができるかどうかを試してみたいと思います。

参考
Getting started — PyStan 2.14.0.0 documentation

はじパタ第1章はじめに

　平井有三先生の「はじめてのパターン認識」を読み始めた。パッと目次を見た感じだとだいたい知ってることが多そうだったが、基礎は大事なので丁寧に読み進めたい。
　今日は第1章「はじめに」を読んだ。100次元の超立方体は毬栗のような形をしているという話がおもしろかった。
　以下に理解度チェックリストを作っておく。

よいパターン認識とは何かいかのキーワードを使って述べよ。

特徴抽出
識別規則の学習
般化能力

以下の尺度の定義、具体例を述べよ。

名義尺度
順序尺度
間隔尺度
比例尺度

ダミー変数とはなにか？
次元の呪いとは？
d次元超立方体について以下を計算せよ。

頂点の数はいくつ？
辺の数は？
超立方体を構成するm次元超平面(1 < m < d)の個数はいくつ？

2017年3月29日水曜日

緑本第11章空間構造のある階層ベイズモデル

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。ようやく最終章までたどり着いた！第11章「空間構造のある階層ベイズモデル」を読んだので、理解度チェックリストを作っておく。

空間相関とは？

どのようなときに考慮しないといけない？
どのように表現できる？

条件つき自己回帰モデル（conditional auto regressive, CAR）とは？
確率場とは？

マルコフ確率場とは？
ガウス確率場とは？

観測データに欠測がある場合、空間相関を考慮する・しないで予測結果はどのように変わる？

2017年3月26日日曜日

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第10章「階層ベイズモデル -GLMMのベイズモデル化-」を読んだ。

　そもそもこの本を読み始めたのは、MCMCをどのようにベイズに適用するのかを理解するのが目的だったのだが、この章に来てようやくイメージができた気がする（9章からぼんやりわかり始めていたものが10章でしっかりと分かった）。この本を読む前までは、MCMCというのを最適化の手法のように勘違い（おそらくメトロポリス法が最適化手法の一つであるsimulated annealingに似ていたのでそう思ってしまったのかも）していた気がする。つまり最尤推定でやっているようにMCMCで対数尤度関数を最適化するパラメータを探すのだと勘違いしていて、それってどうやるのか全然わからない、っていうかわざわざそんなことする意味あるの？と悩んでいたのである。
　しかし、MCMCを「ベイズ」に適用するので、やろうとしていることはパラメータの点推定ではなく、パラメータの分布推定である。パラメータの分布が分かれば、特定パラメータの周辺分布や、応答変数の信頼区間や、ある説明変数が応答変数に影響するかどうかなどの統計検定ができる。こういうことをやるためにMCMCを使うというのが分かると、あぁなるほどな〜と使い方（ベイズのモデルをMCMCで計算するという意味）をイメージできるようになった。

　まえおきが長くなってしまったが、いつもの理解度をチェックするための確認リストを作っておく。

階層ベイズモデルとは？

GLMMだと何が難しい？

階層事前分布とは？
超事前分布とは？
ベイズ統計モデルで主観的な事前分布の使用を避けたい場合は、どのような分布を使えばよいか？
無情報事前分布と階層事前分布のどちらを使うべきかはどのように判断すればよいか？

大域的なパラメータと局所的なパラメータ

緑本第9章 GLMのベイズモデル化と事後分布の推定

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。第9章「GLMのベイズモデル化と事後分布の推定」を読んだので理解度チェックリストを作っておく。

無情報事前分布とは？

一様分布以外にはどのようなものが使える？

データの中央化とは？
データの標準化とは？
MCMCでサンプリングする際に設定すべき値は？

どのようにチューニングするのがよいか？

周辺事後分布とは？
ギブスサンプリングのアルゴリズムについて説明してください

全条件つき分布とは？
メトロポリス法と比較したときの利点は？

共役事前分布とは？

2017年3月25日土曜日

緑本第8章マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。8章「マルコフ連鎖モンテカルロ（MCMC）法とベイズ統計モデル」を読んだので、自分の理解度確認のためのチェックリストを作っておく。

統計モデルのパラメータqの最尤推定値を求めたいとします。qを解析的に求められないときにはどのような方法を使えばよいか？原始的な方法を説明してください。
メトロポリス法のアルゴリズムを説明してください。
MCMCとは？

どの部分がマルコフ連鎖なのか？
どの部分がモンテカルロ法なのか？
MCMCの目的は？

何を出力する？
どういうときに重宝する？

MCMCの定常分布は何を表す？
頻度主義とベイズ統計学の違いについて説明せよ。

モデリングにデータを当てはめた結果えられるものはそれぞれ何か？

2017年3月19日日曜日

緑本第7章一般化線形混合モデル -個体差のモデリング-

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は、第7章「一般化線形混合モデル -個体差のモデリング-」を読んだ。7章からは本の後半となり、より応用的な内容になる。
　以下は自分の理解度をチェックするために作った確認リスト。

過分散とは？
過分散が起こる原因は？
一般化線形混合モデル（GLMM）とは？

どのようなときに使うモデル？

一般化線形混合モデルの最尤推定はどのように行う？

個体差の分布はどのように表すか？
混合分布とは？

データの取り方において、反復、疑似反復の違いは？

GLMMが必要なのはどっち？

2017年3月18日土曜日

緑本第6章 GLMの応用範囲をひろげる -ロジスティック回帰など-

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第6章「GLMの応用範囲をひろげる -ロジスティック回帰など-」を読んだ。
　以下は自分の理解度をチェックするために作った確認リスト。

以下のようなデータをGLMで分析したい場合に、確率分布は何を用いるのがよさそうか？

上限のないカウントデータ
上限のあるカウントデータ
確率変数のとりうる範囲が0以上の連続値データ
負の値もとりうる場合の連続値データ

以下の文のX, Y, Zに当てはまる語句は何か？

ロジスティック回帰は確率分布関数に(X)を、リンク関数に(Y)関数を用いたGLMである。
(Y)関数は(Z)関数の逆関数である。

オッズ比とは何か？
ロジスティック回帰モデル（相互作用項無し）において、ある説明変数xが1増加した場合、オッズ比はどのように変化するか？
交互作用項とは何か？
交互作用項をいれる場合の注意点は？
観測データの割算値を使って回帰を行う場合の問題点は？
オフセット項とは何か？そのメリットは？

2017年3月14日火曜日

緑本第5章 GLMの尤度比検定と検定の非対称性

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第5章「GLMの尤度比検定と検定の非対称性」を読んだ。
　以下は自分の理解度をチェックするために作った確認リスト。

帰無仮説とは？
対立仮説とは？
有意水準とは？
Neyman-Pearsonの検定の枠組みとは？
尤度比検定について説明してください。
第1種の過誤と第2種の過誤の違いは？

Neyman-Pearsonの検定の枠組みではどちらに専念する？

P値とは？
以下の尤度比検定を行うための手法について説明してください。

パラメトリックブートストラップ法
カイ二乗分布を使った近似

「帰無仮説が棄却できなかったので、帰無仮説は正しい」という主張を否定してください。

Neyman-Pearsonの検定の枠組みの非対称性

検定力とは？
検定とモデル選択の目的の違いを説明してください。

2017年3月12日日曜日

緑本第4章 GLMのモデル選択 -AICとモデルの予測の良さ-

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第4章「GLMのモデル選択 -AICとモデルの予測の良さ-」を読んだ。

　以下は自分の理解度をチェックするために作った確認リスト。

モデル選択において、「線形モデルの場合は決定係数R²は高ければ高いほどよい。」という主張を否定してください。
逸脱度とは何でしょうか？

最大対数尤度を使って式で書いてみてください。

フルモデル、ヌルモデルについて説明してください。

それぞれのパラメータ数
最小逸脱度、最大逸脱度

残差逸脱度とは何でしょうか？
AICとは何でしょうか？

逸脱度Dとパラメータ数kを使って式で書くとどうなる？

平均対数尤度とは何でしょうか？
AICでモデル選択をしていい理由を以下の語句を使って説明してください。

最大対数尤度
平均対数尤度
バイアス補正

2017年3月11日土曜日

緑本第3章一般化線形モデル（GLM）- ポアソン回帰 -

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第3章「一般化線形モデル（GLM）- ポアソン回帰 - 」を読んだ。

　以下は自分の理解度をチェックするために作った確認リスト。

ポアソン回帰について説明せよ。

どのようなデータに対して適用するのが自然か？
モデルを式で書くとどうなる？

一般化線形モデルを用いる際に指定する3つの項目を答えよ。
GLMにおいて因子型の説明変数はどのように扱えばよいか説明せよ。
対数リンク関数と恒等リンク関数の違いについて説明せよ。

各説明変数が応答変数に与える影響はどうなるか？

直線回帰と一般化線形モデルの関係について説明せよ。
直線回帰による限界を具体例を用いて説明せよ。

直線回帰が前提としている2つの仮定は何か？

2017年3月10日金曜日

緑本第2章確率分布と統計モデルの最尤推定

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読んでいる。今日は第2章「確率分布と統計モデルの最尤推定」を読んだ。

　以下は自分の理解度をチェックするために作った確認リスト。

確率分布とは何か？
以下の確率分布はどのような場面で使うのが適切か？

ポワソン分布
二項分布
正規分布
ガンマ分布

尤度とは何か？
以下のデータはポワソン分布にしたがってランダムに抽出されたものとする。

4, 3, 4, 2, 2, 2, 1, 3, 6, 3
対数尤度を最大化するパラメータと標本平均が一致することを確認せよ。（ポワソン分布の式をそらで書けないようだとダメ）

標準誤差とは何か？
推定と予測の違いについて説明せよ。

2017年3月9日木曜日

緑本第1章データを理解するために統計モデルを作る

　久保先生の緑本こと「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC」を読み出した。自分自身の理解度チェックのために各章ごとに確認リスト的なものを作ることにした。確認リストが質問形式なのは何となく詰め寄られてる気がして緊張感が出るから。

　今日は第1章「データを理解するために統計モデルを作る」を読んだ。

統計モデルとは何でしょうか？
応答変数と説明変数について説明してください。
リンク関数について説明してください。
統計モデルのどのあたりが「統計」なモデルなのか説明してください。
統計学の誤用あるあるネタを挙げてください。

登録: 投稿 (Atom)

Tech Tips

Page List

Search on the blog