ベイズ統計学

HOME
１．予測理論
ベイズ統計学
Stanによる統計モデル：推測統計の基本から時系列分析まで

2017年12月2日 / 最終更新日時 : 2017年12月2日馬場真哉ベイズ統計学

Stanによる統計モデル：推測統計の基本から時系列分析まで

新規作成：2017年12月2日
最終更新：2017年12月2日

StanとRを用いた統計モデル構築の基本について説明します。
統計学の初歩からベイズ推論、Stanというソフトウェアの概要といった基本事項から時系列モデルの推定の方法まで、順を追って説明します。

この記事はStan Advent Calendar 2017の2日目の記事となります。
詳細は「Stan Advent Calendar 2017」を参照してください。このリンクをたどると、Stanに関する様々な記事を読むことができます。

スポンサードリンク

1．推測統計学の基本

たとえば「３」というデータが私たちの手元にあったとします。
ここで「あぁ、３があるなぁー」と思うのが推測統計学を知らない人です。
ここで「３というデータはどのようなプロセスで私たちの手元にやってきたのだろう」と考えるのが推測統計学です。

３というデータが、例えば魚の体長のデータだったとしましょう。近所の小さな池で釣りをしたら、３センチの魚が釣れたんだそうです。
ここで「あぁ、３センチの魚が釣れたんだなぁー」と思うのが推測統計学を知らない人です。
ここで「３センチの魚が釣れたということは、湖の中に３センチの魚がいたということだ」と考えるのが推測の第一歩です。
そのうえで「湖の中には３センチの魚がいた。次に釣りをしたらやっぱり同じくらいの大きさの魚が釣れるんじゃないだろうか」と推測するのが推測統計学というものです。

推測統計学では、データが得られるプロセスを抽象化します。
ここで登場するのが確率分布です。

近所の小さな池の中に、10尾の魚が住んでいたとしましょう。
10尾の魚の大きさはすべてわかっていて、以下の通りだったとします。
1センチ：1尾
2センチ：2尾
3センチ：4尾
4センチ：2尾
5センチ：1尾

この時、池から魚を1尾捕まえて「体長というデータが得られる」というプロセスは、以下の確率分布から乱数を1つ生成するプロセスと同じだとみなすことができます。
1センチ：10%
2センチ：20%
3センチ：40%
4センチ：20%
5センチ：10%

ほかにも、例えば、赤い魚が1万尾、青い魚が1万尾の合計2万尾の魚が住んでいるちょっと大きな湖から「データが得られるプロセス」を考えてみましょう。
これは以下の確率分布から乱数を1つ生成するのと同じです。
赤い魚：50%
青い魚：50%

これは「コインを投げて表が出たら赤い魚、裏が出たら青い魚」だとみなすことで簡単に「データが得られるプロセス」をシミュレートすることができます。

こういう確率分布を「母集団分布」と呼びます。
推測統計学では、データが得られるプロセスとして確率分布を使うということを是非覚えておいてください。

推測統計学では、この母集団分布を頑張って推定します。
母集団分布が推定できれば「次に魚を釣ったらどんな魚が釣れるだろうか」を予測することもできるようになりますね。

ところで、母集団分布を、何もわからない状況から推定するのはちょっと大変です。
そこで正規分布などの「特徴が良く知られた確率分布」を使うことが多いです。
もしも母集団分布に正規分布を仮定したら、平均と分散という2つのパラメタ(母数と呼ばれます)を推定するだけで母集団分布の推定が完了します。

別に正規分布じゃなくてもいいですし、正規分布を使うと現実の「データが得られるプロセス」と乖離してしまうこともあります。こういう時はポアソン分布など別の確率分布を使う必要があります。
この記事では、まずは簡単のため、母集団分布として正規分布を仮定して進めていきます。

2．ベイズ推論の基本

母集団分布が推定できると嬉しいな、という話をしました。
次は母集団分布を推定する方法論を説明します。
データと「ベイズの定理」を使って母集団分布を推定します。

ベイズの定理あるいはベイズ更新の考え方については「ベイズ統計学基礎」を参照してください。この記事では計算の方法に絞って解説します。
また、確率密度という考え方も知っておくと、この記事が読みやすいかと思います。詳しくは「確率密度関数と正規分布」を参照してください。この記事では確率密度と確率を特に区別しないこともあります。

問題を整理します。
母集団分布としては正規分布を仮定しています。
そのため、あとは平均と分散という2つのパラメタを推定することで、母集団分布を推定することができます。

両方をいっぺんに扱うのは少し難しいので、まずは平均値の推定を考えます。
これは「母集団の平均値」すなわち母平均を推定する問題だと考えてもらっても結構です。しかし、後々統計モデルを学ぶことになるので「正規分布のパラメタ(母数)」を推定しているのだという認識を持ってもらえると嬉しいです。

データがない状況で母平均はいくらになるかと聞かれても困るんじゃないかと思います。
３かもしれないし６３７かもしれない。この「よくわからない」という状況を表すのに便利なものが無情報事前分布です。
無情報事前分布とは、例えば「平均0、分散1000000」の正規分布などが使われます。
母平均が無情報事前分布に従うと考えると、「母平均が3になる確率密度」も「母平均が637になる確率密度」も共に(あるいは平等に)とても小さな値になります。
３になるのか６３７になるのか－３５になるのか見当もつかないというときは、みな等しく平等に低い確率を割り当てた確率分布を使うということです。

ここで「３」というデータが手に入ったとしましょう。
そうしたらベイズの定理を使って無情報事前分布を更新します。「３」というデータが手に入ったのだとしたら、平均値は３になりやすくなると考えられますね。そのため「母平均が３になる確率」がちょっと増えることになります。逆に「母平均が６３７である確率」や「母平均が－３５である確率」は減ります。
更新された後の母平均の確率分布は、事後分布とも呼ばれます。
標語のように言うと、ベイズの定理は事前分布を事後分布に更新する数式だと思えばよいです。

まとめます。
データは母集団分布から得られたと考えます。
母集団分布に正規分布を仮定すると、正規分布を構成するパラメタ、すなわち母平均と母分散が推定できれば、母集団分布が推定できることになります。
例えば母平均を推定することを考えます。
データがなければ、母平均がいくつになるのかさっぱりわかりません。そこで母平均を無情報事前分布に従う変数だと考えます。そうすることで「ありうる母平均の値」にみな等しく平等に低い確率を割り当てることができます。こうしておけば「母平均がいくらなのかよくわからない」という状況を数学的に表現できます。
データが手に入ったらベイズの定理を使って事前分布を事後分布に更新します。例えば「３」というデータが手に入ったならば「平均値は３でありやすい」と更新されることになります。
母平均が「３でありやすい」ことがわかれば、私たちが手に入れたそのデータは「母平均が３の正規分布から得られたのではないか」と推測することができるようになります。そして「次に得られるデータも、３に近い値かもしれない」と予測することができるでしょう。

ここで少し注意点を述べます。
2つの確率分布が出てきているので、それをごっちゃにしないように気を付けてください。
データが得られるプロセスとしての母集団分布がその1つです。
もう1つは「母集団分布のパラメタがいくらになるのかよくわからない」という「わからなさ」を定量化するために使われている確率分布です。

母集団分布のパラメタがいくらになるのかよくわからないという「わからなさ」を定量化するために「母集団分布のパラメタが従う確率分布」を想定しているということです。

3．StanとMCMC法の基本

MCMC法とは、任意の確率分布に従う乱数を生成する手法です。
ギブスサンプラーやハミルトニアンモンテカルロなどいろいろなアルゴリズムが提案されています。

Stanとはハミルトニアンモンテカルロ法を用いた乱数生成の機能を持つソフトウェアです。

MCMC法を使うことで、事後分布に従う乱数を生成することができます。
事後分布に従う乱数のヒストグラムを描けば、それが事後分布の形状となりますね。
事後分布の平均値を求めたいと思ったら、MCMC法により生成された「事後分布に従う乱数の平均値」を計算すればよいです。
事後分布の確率密度関数はとても複雑なものになることがしばしばあります。そんな時でもMCMC法を使えば、事後分布に従う乱数を相手にするだけで、推定結果の解釈ができるので、簡単です。

4．Stanによる平均と分散の推定

Stanを使って母集団分布を推定してみましょう。
やることは以下の3つです。
①母集団分布の構造を決める(今回は正規分布を仮定するので、平均と分散という2つのパラメタを推定すればOKです)
②MCMC法を使って、母集団のパラメタの事後分布に従う乱数を生成する
③生成された乱数を使って、推定結果を解釈する

Stanの細かいコードの説明などは「Stanによるベイズ推定の基礎」も併せて参照してください。ソフトのインストールの方法などもこちらに書いてあります。

まずは分析のためのライブラリを読み込みます(Stanのインストールはすでに終わっているとしています)。

# ライブラリの読み込み
library(rstan)

# 計算の並列化
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

データを用意します。
池で釣れた魚の大きさという架空のデータを使います。
『fish_data.csv』

データを読み込みます。

> fish_data <- read.csv("https://logics-of-blue.com/wp-content/uploads/2017/12/fish_data.csv")
> head(fish_data)
         x
1 2.373546
2 3.183643
3 2.164371
4 4.595281
5 3.329508
6 2.179532

ここからベイズモデリング開始ですね。
①母集団分布の構造を決める(今回は正規分布を仮定するので、平均と分散という2つのパラメタを推定する)作業に移ります。

母集団分布の構造は、Stanファイルに記述します。
以下のコードを『estimate-mean.stan』という名称で保存します。
Stanファイルは作業ディレクトリに保存するように注意してください。
作業ディレクトリは、Rで「getwd()」と書いて実行すればわかります。

data {
  int N;
  real length_data[N];
}

parameters {
  real mu;               // 母平均
  real<lower=0> sigma;   // 母分散
}

model {
  for(i in 1:N)
    length_data[i] ~ normal(mu, sqrt(sigma));
}

Stanファイルには最低でもdata、parameters、modelの3つのブロックが必要となります。

dataブロックには、名前の通り分析対象となるデータの構造を指定します。
ここではサンプルサイズNと魚の体長データlength_dataを指定しました。
『int N』とすることでサンプルサイズNは整数として扱われます。
『real length_data[N]』とすることで、length_dataは長さがNの実数として扱われます。

parametersブロックには、推定すべきパラメタを指定します。
今回は母平均と母分散を推定するため、この2つを指定しました。
<lower=0>とすることで、最小値が0であると指定することができます。

modelブロックには「データが得られるプロセス」を直接指定します。
今回はN個の体長データが「平均mu、標準偏差sqrt(sigma)の正規分布」に従って得られると仮定しているので、このように書きます。
「データが○○という確率分布に従う」というのは「データ　～　確率分布」と表記することは覚えておくと良いと思います。
(ベクトル化をするともう少し短くコードを書くことができますが、コードが見にくくなるので、この記事では使いません。)

モデルの構造が決まりました。
「②MCMC法を使って、母集団のパラメタの事後分布に従う乱数を生成する」の作業に移ります。

ここからはまたRに戻ります。
まずはStanに渡すデータを作ります。これはStanファイルのdataブロックと対応している必要があります。

> # Stanに渡すデータ
> N <- 200
> data_stan_1 <- list(length_data = fish_data$x, N = N)
> data_stan_1
$length_data
  [1] 2.3735462 3.1836433 2.1643714 4.5952808 3.3295078 2.1795316 3.4874291 3.7383247 3.5757814
 [10] 2.6946116 4.5117812 3.3898432 2.3787594 0.7853001 4.1249309 2.9550664 2.9838097 3.9438362
 [19] 3.8212212 3.5939013 3.9189774 3.7821363 3.0745650 1.0106483 3.6198257 2.9438713 2.8442045
 [28] 1.5292476 2.5218499 3.4179416 4.3586796 2.8972123 3.3876716 2.9461950 1.6229404 2.5850054
 [37] 2.6057100 2.9406866 4.1000254 3.7631757 2.8354764 2.7466383 3.6969634 3.5566632 2.3112443
 [46] 2.2925048 3.3645820 3.7685329 2.8876538 3.8811077 3.3981059 2.3879736 3.3411197 1.8706369
 [55] 4.4330237 4.9803999 2.6327785 1.9558654 3.5697196 2.8649454 5.4016178 2.9607600 3.6897394
 [64] 3.0280022 2.2567268 3.1887923 1.1950414 4.4655549 3.1532533 5.1726117 3.4755095 2.2900536
 [73] 3.6107264 2.0659024 1.7463666 3.2914462 2.5567081 3.0011054 3.0743413 2.4104791 2.4313313
 [82] 2.8648214 4.1780870 1.4764332 3.5939462 3.3329504 4.0630998 2.6958161 3.3700188 3.2670988
 [91] 2.4574800 4.2078678 4.1604026 3.7002136 4.5868335 3.5584864 1.7234078 2.4267346 1.7753874
[100] 2.5265994 2.3796333 3.0421159 2.0890784 3.1580288 2.3454154 4.7672873 3.7167075 3.9101742
[109] 3.3841854 4.6821761 2.3642635 2.5383553 4.4322822 2.3493036 2.7926193 2.6071921 2.6800071
[118] 2.7208867 3.4941883 2.8226695 2.4940425 4.3430388 2.7854206 2.8204435 2.8998093 3.7126663
[127] 2.9264356 2.9623658 2.3183395 2.6757297 3.0601604 2.4111055 3.5314962 1.4816059 3.3065579
[136] 1.4635502 2.6990239 2.4717201 2.3479052 2.9431032 1.0856406 4.1765833 1.3350276 2.5364696
[145] 1.8840799 2.2491810 5.0871665 3.0173956 1.7136995 1.3593945 3.4501871 2.9814402 2.6819316
[154] 2.0706379 1.5125397 1.9248077 4.0000288 2.3787333 1.6155732 4.8692906 3.4251004 2.7613529
[163] 4.0584830 3.8864227 2.3807570 5.2061025 2.7449730 1.5755053 2.8556004 3.2075383 5.3079784
[172] 3.1058024 3.4569988 2.9228471 2.6659992 2.9652740 3.7876396 5.0752450 4.0273924 4.2079084
[181] 1.7686766 3.9838956 3.2199248 1.5327500 3.5210227 2.8412454 4.4645873 2.2339180 2.5697882
[190] 2.0738905 2.8228960 3.4020118 2.2682518 3.8303732 1.7919172 1.9520156 4.4411577 1.9841525
[199] 3.4119747 2.6189239

$N
[1] 200

次にMCMCを実行します。

fit_stan_1 <- stan(
  file = "estimate-mean.stan",
  data = data_stan_1,
  seed = 1
)

stan関数を使って「パラメタの事後分布に従う乱数」を生成します。
stan関数には、Stanファイルの名称、Stanに渡すデータ、そして乱数の種を指定します。
乱数は文字通り「ランダム」に発生しますが、乱数の種を指定すると、毎回同じ値が出るようになります。分析の再現性を担保するために、seedはなるべく指定します。
計算が終わるのにはちょっと時間がかかります(1分ほどです)。

「③生成された乱数を使って、推定結果を解釈する」に移ります。
結果はprint関数で確認します。

> # 推定結果の表示
> print(fit_stan_1, digits = 2)
Inference for Stan model: estimate-mean.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
mu      3.04    0.00 0.07   2.91   2.99   3.03   3.08   3.16  3072    1
sigma   0.88    0.00 0.09   0.73   0.82   0.87   0.94   1.07  2959    1
lp__  -85.91    0.02 0.96 -88.46 -86.31 -85.63 -85.22 -84.96  1788    1

Samples were drawn using NUTS(diag_e) at Wed Nov 29 15:58:22 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).

これを見ると、MCMCにより作成された「事後分布に従う乱数」の平均値や95%区間がわかります。
右端の『Rhat』は、乱数が正しく生成できたかどうかをチェックする指標です。この値が1.1未満であればOKです。

母平均muは95%の確率で2.91から3.16の間にあるだろうと推定されました。
この幅が広ければ「母平均の値がいくらなのかよくわからない」ということになります。逆に狭い範囲であれば母平均の値に関して信頼できる推定値が出せたということになります。
母平均の値を1つ推定値として提示しろ、といわれれば、生成された乱数の50%点(中央値)である『3.03』と答えます(生成された乱数の平均値である『3.04』でも構いませんが、中央値の方が多く使われる印象があります)。

何はともあれ、これで母集団分布の推定が完了ですね。
母集団分布は「平均2.99、分散0.87の正規分布」であると推定できたわけです。
次に同じ池から同じ条件で釣りをしたら、やっぱり「平均2.99、分散0.87の正規分布」からデータが得られると想定できます。

スポンサードリンク

5．データ生成過程

応用編として、時系列モデルを構築してみます。
時系列分析については「時系列解析_理論編」なども参照してください。

時系列モデルを構築する大きな目的は、データ生成過程を推定することです。
時間によって変化する確率分布をデータ生成過程と呼びます。

この記事では自己回帰モデルと呼ばれるデータ生成過程を対象とします。
1次の自己回帰モデルは以下のように表記されます。
1次というのは「1時点前までの情報を使った」くらいの意味です。

データ～ β × 1時点前のデータ＋ノイズ

数式で書くと以下のようになります。
$$
y[i] \sim gauss(\beta \cdot y[i – 1], \sigma) \\
$$

βは回帰係数です。添え字の[i]は時点を表すインデックスです。
gauss()は、正規分布です。平均が「β × 1時点前のデータ」であり分散がσの正規分布に従ってデータが得られると仮定しています。

6．Stanによる自己回帰モデル

stanを使って自己回帰モデルを推定してみます。
stanファイルは以下のようになります。「autoregressive.stan」という名前で保存します。


data {
  int N;
  real y[N];
}

parameters {
  real beta;
  real<lower=0> sigma;
}

model {
  for(i in 2:N)
    y[i] ~ normal(beta * y[i-1], sqrt(sigma));
}

自己回帰モデルでは、回帰係数betaと、分散の大きさsigmaを推定します。
i時点のデータは、『beta ×「i-1時点のデータ」＋ノイズ』として得られると考えています。

分析対象となるデータはシミュレーションをして作ります。
真のデータ生成過程は以下の通りです。

y ～ 0.5 × 1時点前のy＋分散1の正規分布に従うノイズ

> # シミュレーションデータの作成
> set.seed(1)
> sim_data <- arima.sim(n = N, model = list(ar = c(0.5)))
> 
> # stanに渡すデータ
> data_stan_2 <- list(y = sim_data, N = N)
> data_stan_2
$y
Time Series:
Start = 1 
End = 200 
Frequency = 1 
  [1]  1.614242003  1.196964238 -0.022758462 -2.226079118  0.011891359 -0.038987929 -0.035684228
  [8]  0.925994097  1.284218243  1.236010443  1.536982593  1.550627597  0.849878782 -1.564412305
 [15] -0.162380405 -0.137318942 -0.224454978 -1.582979873 -1.269639991 -0.216878436  1.250240334
 [22]  0.522332440  0.648837831  0.270613875 -1.241752619 -1.035870873 -0.912225390 -0.515426092
 [29]  0.842312326  1.184331911  0.427642359 -0.039540500  0.677193125  0.895259761 -0.241125814
 [36] -0.828058064 -0.049447070  0.743809390  0.259558483  1.010886968  0.903549364 -0.160251711
 [43]  0.260993836 -0.998866178  0.933590613  2.447195205  0.856376126 -0.615946563  0.261746346
 [50] -0.004181431  2.399527045  1.160523520  1.270001122  0.663002720 -0.411771849 -0.017093625
 [57] -1.813505441  0.558802141  0.432654409  2.388938875  1.669978966  0.125043052  0.673247880
 [64] -0.597473692 -1.552370246 -0.484738888 -0.685661317 -0.341725307 -0.096521329 -0.637781611
 [71] -0.887559538 -0.578958384  0.888607804 -1.079262898  0.054314739  0.360107740  1.243153708
 [78]  0.317392930  0.528715275  0.531456428 -0.276791817  1.069471898  1.695138564  1.547782932
 [85]  2.360724920  1.738848886 -0.407167766 -0.776849297 -1.613037263 -1.279919268 -1.260326311
 [92] -0.588047283 -1.204945290 -0.444443872 -0.876806580  1.328883979  1.381149466  1.600748962
 [99]  1.184559839  2.274456000  0.501491546 -0.210898957  1.326832760  0.012720027 -0.201020730
[106] -0.493318295 -0.566652016 -0.562439311  0.212968676 -0.070846144 -0.541380534  1.072348558
[113]  0.321594870 -0.018759095 -0.109570289  0.657881163  0.255376177  0.090053917 -0.636633520
[120] -0.642587032 -0.261133076 -0.719461024  0.171765681 -1.432511242 -0.409697760 -1.741298704
[127] -1.171625479 -1.114092644 -1.209141103 -0.661467329 -2.245093090  0.054036767 -1.637954053
[134] -1.282507428 -1.757173819 -1.629405911  1.272463590  0.653627415 -0.959486823 -2.120348946
[141] -0.609987372 -0.323553519 -0.479845134 -1.169284714 -2.072102667 -2.111243630 -0.055593011
[148] -0.649063201 -1.708958448  1.014811399  0.932506077  0.227605937  1.172286017  1.472565660
[155]  0.117039782  2.264622355  0.877284148 -0.985852576 -0.637325890 -0.111124606  2.252416096
[162]  1.232010416  1.073004013  0.459349071 -0.104326307 -0.086889182  0.744195015  2.447342516
[169]  2.251063697  2.333440247 -0.064603298  0.951593921  0.695721764 -1.119389147 -0.038671831
[176] -0.178090520  1.375542052 -0.078310974 -0.469367241 -1.160793118 -0.757500520  0.023261519
[183] -0.720117413  0.470314461 -0.972925556 -1.534447191  0.673934112 -0.678880410  0.072534508
[190] -0.344808797  0.236997441  1.807372007  2.490274437  0.914229418 -1.828120826  1.583601177
[197]  1.458866755  1.270760714  0.621980834  0.821098840

$N
[1] 200

MCMCを実行します。

# MCMCの実行
fit_stan_2 <- stan(
  file = "autoregressive.stan",
  data = data_stan_2,
  seed = 1
)

結果はこちら。余計な出力は省略します。

> # 推定結果の表示
> print(fit_stan_2, digits = 2)
Inference for Stan model: autoregressive.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

        mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
beta    0.48    0.00 0.06   0.36   0.44   0.48   0.52   0.60  3675    1
sigma   0.94    0.00 0.10   0.77   0.87   0.93   1.00   1.14  3346    1
lp__  -91.74    0.02 0.97 -94.21 -92.15 -91.45 -91.05 -90.77  1768    1

回帰係数betaは、95%の確率で0.36から0.60の間に収まるということになりました。
分散sigmaは95%の確率で0.77から1.14の間に入ると解釈できます。

7．Stanによる一般化自己回帰モデル

先ほどは母集団分布として正規分布を仮定していました。
これは、－∞～＋∞までの値をとる連続的な変数に対してしか用いることができません。

例えば、商品が売れた個数や観測された毎日の小鳥の数などは0以上の整数しかとりません。
こういった場合は、母集団分布に、正規分布ではなくポアソン分布を使います。
ポアソン分布の母数は平均値λ(ラムダ)のみです。

λは負の値になってはいけません。
そのため、指数関数を途中にかませます。

シミュレーションをしてデータを作ります。

> # シミュレーションデータの作成
> set.seed(1)
> sim_pois <- numeric(N)
> for(i in 1:N){
+   sim_pois[i] <- rpois(n = 1, lambda = exp(sim_data[i]))
+ }
> 
> # stanに渡すデータ
> data_stan_3 <- list(y = sim_pois, N = N)
> data_stan_3
$y
  [1]  4  3  1  1  0  2  3  3  4  1  3  3  3  0  1  1  1  2  0  1  6  1  2  0  0  0  0  0  4  2
 [31]  1  1  2  1  2  1  2  0  2  2  4  1  2  0  2 14  2  1  2  1 14  2  1  0  0  1  0  1  3  9
 [61]  4  1  1  0  0  0  1  0  2  0  0  0  4  1  1  2  7  1  2  1  0  4  3  6  6  4  0  1  1  1
 [91]  1  0  0  1  0  4  3  3  8 11  1  0  3  3  1  3  1  0  1  0  0  4  0  1  1  6  1  1  0  1
[121]  0  0  0  0  1  0  0  0  0  2  0  1  0  2  0  0  4  1  0  0  0  0  1  0  0  0  1  0  0  3
[151]  1  1  1  3  0  8  4  0  1  2  9  1  2  2  0  1  4 15  8 14  1  3  4  0  0  2  4  2  0  1
[181]  1  3  0  2  0  0  4  0  1  0  2  5 14  1  0  3  2  3  2  1

$N
[1] 200

自己回帰モデルをポアソン分布に従うデータに適用させます。
以下のようにモデル化することになります。
$$
\mu[i] \sim gauss(\beta \cdot \mu[i – 1], \sigma) \\
\lambda[i] = \exp ( \mu[i] ) \\
y[i] \sim pois(\lambda[i])
$$

この時のStanファイルは以下のようになります。「autoregressive_pois.stan」という名称で保存します。

data {
  int N;
  int y[N];
}

parameters {
  real beta;
  real<lower=0> sigma;
  real mu[N];
}

transformed parameters{
  real lambda[N];
  
  for(i in 1:N)
    lambda[i] = exp(mu[i]);
}

model {
  for(i in 2:N)
    mu[i] ~ normal(beta * mu[i-1], sqrt(sigma));
    
  for(i in 1:N)
    y[i] ~ poisson(lambda[i]);
}

変わったのは2点。
1つはtransformed parametersブロックが追加されました。
平均値λが0以上な必要があるため、指数関数exp()をかませるために使いました。

もう1つはmodelブロックです。
データ(y)は平均がλのポアソン分布に従うように変更しました。

MCMCを実行します。

# MCMCの実行
fit_stan_3 <- stan(
  file = "autoregressive_pois.stan",
  data = data_stan_3,
  seed = 1
)

推定結果はこちら。
結果が長いので『pars = c(“beta”, “sigma”)』として表示するパラメタを指定しました。

> # 推定結果の表示
> print(fit_stan_3, digits = 2, pars = c("beta", "sigma"))
Inference for Stan model: autoregressive_pois.
4 chains, each with iter=2000; warmup=1000; thin=1; 
post-warmup draws per chain=1000, total post-warmup draws=4000.

      mean se_mean   sd 2.5%  25%  50%  75% 97.5% n_eff Rhat
beta  0.47    0.00 0.11 0.25 0.39 0.47 0.54  0.66  1326    1
sigma 0.90    0.01 0.18 0.59 0.77 0.89 1.01  1.30   973    1

betaがおよそ0.5でsigmaがおよそ1前後となっているので、正しく推定されていることがわかります。

8．モデルを組む時に考えていること

stanによるモデリングは、Rのglm関数やar関数を使った時などと比べると、様々な”改造”ができるため、自分の思い通りのモデルを構築しやすいと言えます。
しかし、書くべきコードの量が増えますし、そもそも、どういうモデルを想定すればいいのかをちゃんと自分で考えてからでなければモデル化ができません。
自由なモデリングっていうのは、裏を返せば、どういうモデルにするかを全部ユーザー任せにしたということですので。

ここではStanのコードをどうやって「思いつく」のかを、管理人個人の例ではありますが、簡単に言葉で説明してみます。

最も重要なことは、手に入ったそのデータが「何らかの確率分布に従う確率変数」であることを理解することだと思います。
例えば、手持ちの「３」というその数値は「平均3、分散1の正規分布」に従う確率変数だ、とかそういう認識です。
そのうえで「データ～　gauss(3, 1)」のような表記の仕方に慣れることができれば、最初の関門クリアです。

統計モデルを構築する作業を
1．母集団の確率分布を選ぶ(正規分布とかポアソン分布とか)
2．確率分布の母数の変化のパターンを数式で表現する
の2つに分けて考えると話が簡単になります。
自己回帰モデルですとまずは母集団分布に正規分布を仮定しました。
そのうえで「前の時点のデータが大きな値だったら、次の時点も大きなデータとなりやすい」という状況を、「正規分布の平均値が1時点前と関係している」という状況だと考え、「データ～ gauss(β × 1時点前のデータ, σ)」とモデル化しました。
頭で想像したモノから「データ～ gauss(β × 1時点前のデータ, σ)」という表現が出てくれば、Stanによるモデル化ができるようになります。

モデルを改造する(例えばポアソン分布に従う自己回帰モデルを組もうと考える)時でも、できれば以下のようにまとめて数式で書けると、Stanコードの見通しがかなり良くなります。
$$
\mu[i] \sim gauss(\beta \cdot \mu[i – 1], \sigma) \\
\lambda[i] = \exp ( \mu[i] ) \\
y[i] \sim pois(\lambda[i])
$$

ここまで書けちゃえば、「＝」でつながった式をtransformed parametersに「～」でつながった式をmodelブロックに入れてしまえば、Stanコードがほとんど完成です。
dataブロックは考えることがほとんどないですし、parametersブロックはmodelブロックに合うように指定すればそれほど難しくないはずです(極端な話、定義し忘れててもエラーになるのでわかる)。

参考文献

ベイズ統計モデリング: R,JAGS, Stanによるチュートリアル原著第2版

とても大きくて分厚い本です。足の上に落とすと痛い(痛かった)です。
分厚い本ですが、逆に言えばそれだけ記述が丁寧だということです。ベイズ推論の考え方からMCMCのアルゴリズムの説明、そしてデータ分析の具体例まで載っています。
初学者の方でもゆっくりならば読み進められる難易度かなと思います。

StanとRでベイズ統計モデリング (Wonderful R)

Stanの使い方が載っている本としては、間違いなく日本で一番詳しい本です。
統計モデルの考え方についての記載もあります。

書籍以外の参考文献

Stan Advent Calendar 2017
　→ この記事も参加しているAdvent Calendarです。Advent Calendarは12月1日～24日まで毎日誰かがその分野の記事を書くというイベントです。Stan以外にもRやPythonなどいろいろあります。「2017年プログラミング言語カテゴリーのカレンダー」も参照してみてください。

Stanの日本語のリファレンス
　→ Stanのマニュアルです。日本語に翻訳してくださっているので簡単に読めます。自己回帰モデルをはじめとした様々なモデルの実装例が載っています。

スポンサードリンク