2016-09-06

RとSTANとでBayesianの基礎

R Stan rstan ベイズ

ベイズ推定をする
- データがある
- モデルがある
- モデルは、パラメタを持っていて、そのパラメタの値を定めると、データを観察する確率が定義できる
- Stanでは、このパラメタがどういう分布を取っているかを乱数を使って標本分布として返す
- パラメタ以外にも、どういう分布になっているかを知りたいものがあれば、それは生成量として標本分布を作ってくれるので、分布がわかる
Stanを実行するには
- rstanパッケージにはstan()という関数があり、それを使ってStanを実行する
- stan()関数を回すのに必須なのは
  - hoge.stan　というテキストファイルと
  - hoge.stan 　に記載された"data"であるRのリストオブジェクト
習うより、慣れろ、の手順
- 問題設定(モデル)を自分で作ってみる
- モデルに合わせて、標本を発生してみる
- モデルに合わせて知りたいことを決めてみる
- モデル、知りたいことに合わせてhoge.stanファイルを作ってみる
- stan()を回してみる
例１
- 問題設定
  - 家の前を通る車が京都ナンバーである割合はどれくらいか知りたくなった。京都ナンバーである割合は未知変数である。これをpとする。ただし、pは0以上1以下の実数
  - 家の前を通る車をn台観察するとする。ただしnは0以上の整数
  - n台のうち、京都ナンバーである車の台数をmとする。ただし、mは0以上n以下の整数
  - p,n,mには確率モデルが入れられる。二項分布になっているものとする。 $Pr(m|n,p) = \begin{pmatrix}n \\ m \end{pmatrxi} p^m (1-p)^{n-m}$
- 標本を発生してみる

p <- 0.7 # 仮にこの値にしてみる
n <- 10 # 仮に10台観察することにする
n.iter <- 1000 # 10台の観察を、もし1000回実行するとしたら
m <- rbinom(1000,n,p)
hist(m)

- 知りたいこと
  - 今、n=10台の観察をして、m=7台が京都ナンバーだったとする
  - 京都ナンバーの割合pの期待値はいくつか
  - 3台ずつ観察してそれを100回繰り返したときに、3台とも非京都ナンバーになっているのは、100回のうち何回かの分布
- hoge.stanを書く

# 京都ナンバーの割合
# dataはRのリストオブジェクトとして取り込む
data { 
  int<lower=0> n;                              #観察車台数
  int<lower=0> m;                              #京都ナンバー車台数
}
parameters {
  real<lower=0,upper=1>   p;            #京都ナンバーの比率(母集団の比率)
}
transformed parameters { #今回は不要な項目
}
model {
  m ~ binomial(n,p);
}
generated quantities{
  int<lower=0> allthree;
  allthree <- 0;
  for(i in 1:100){
    if(binomial_rng(3,1-p) == 3){allthree <- allthree +1;}
  }
}

- 実行結果

n <- 10
m <- 7
out <- stan(file="kyotonumber.stan",dat=list(n=n,m=m))
# 事後分布はstan()関数の出力outにextract()という関数を適用することで抽出できる
out.extracted <- extracted(out)
hist(out.extracted[[1]])
hist(out.extracted[[2]])

例２
- 例１はちょっと簡単すぎる。Stanを回すまでもない…
- 問題設定を変える
  - k台ずつ、複数回(N回)、ナンバープレートを確認したところ、k台とも非京都ナンバーだったのは、N回中n回だった、という情報を得たとする。ここから、例１と同じことが知りたいとする
  - 知りたいことは変わっていないが、変わっているのは、観察事項

# 京都ナンバーの割合
# dataはRのリストオブジェクトとして取り込む
data { 
  int<lower=0> k;                              #観察車台数
  int<lower=0> N;                              #観察セット数
  int<lower=0> n;                              #すべて非京都ナンバーだったセット数
}
parameters {
  real<lower=0,upper=1>   p;            #京都ナンバーの比率(母集団の比率)
}
transformed parameters { #今回は不要な項目
}

model {
  n ~ binomial(N,(1-p)^k);
}
generated quantities{
  int<lower=0> allthree;
  allthree <- 0;
  for(i in 1:100){
    if(binomial_rng(3,1-p) == 3){allthree <- allthree +1;}
  }
}

out2 <- stan(file="kyotonumber2.stan",data=list(k=k,N=N,n=n)