毎回やる儀式

インターネットに接続しておいてください（このサイトが見えていればOK）。
Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
File>New File>R scriptで新しいコードファイルを作ってください。

バッターのデータ，batter.csvをb.dataというオブジェクトに入れているとします。できてないひとは，次の通りやってみよう。

b.data <- read.csv("batter.csv")

今日使うパッケージ

rstan
shinystan

ベイジアンモデリングって？

「ベイジアンモデリング」は「ベイジアンのモデリング」であり，ここでは「ベイズ的に考えるモデリング」というぐらいの意味で使ってます。

「ベイズ統計学」と比較されるのは「頻度主義統計学」などといい，これまで紹介してきたモデリングはほとんどすべて，頻度主義的なアプローチによっていました。頻度主義者とベイジアンの考え方はどう違うのかを簡単にまとめたのが次の表です。

対象	頻度主義者	ベイジアン
母数	定数	確率変数
データ	確率変数	定数

頻度主義者は，目の前のデータはある確率分布から得られたもの＝泡沫の夢，と考えます。このデータを生み出した本当の世界を見つめたいのです。

ベイジアンは，目の前のデータは真実であると考えます。このデータから夢を見たい＝こういうデータが出てくるということは，どういうことが言えるだろうか，という可能性を考えます。

具体的には，頻度主義者にとって回帰係数，パス係数，因子負荷量などモデルの結果として得られたものは「定数」（一つの値）ですので，それを報告するだけです。ベイジアンにとって，回帰係数，パス係数，因子負荷量などモデルの結果として得たいものは「確率変数」ですので，確率分布の平均や分散など，分布の特徴を複数報告する必要があります。

「報告するデータの数が増えるんならベイズは面倒じゃないか」と思うかもしれませんが，ベイジアン・アプローチはいろんな利点があります。例えば

変な仮定(帰無仮説，分散の等質性など)を考えなくても，イメージをそのまま表現できる。ルールに従っているだけじゃなくて，ルールに沿った分布を生み出す側だから。
結構複雑なモデルでも解ける。分布を見つけるんじゃなくて生み出す側だから。
少ないサンプルでも解ける。多くのデータがないと分布としての特徴を当てはめられない，という縛りはなく，多いに越したことはないけど手元のデータから分布を生み出す側だから。

などなど。

ベイジアンモデリングをするには

少数のサンプルしかない場合や，心理的メカニズムをいろいろ組み込んだモデルを考えたい場合，階層的データの上位階層の特徴を表現したい場合など，心理学の多くのシーンでベイジアンモデリングは有効です。ではどうして最初からそれをやらなかったのでしょうか？答えはいつもの通り，時代が追いついてなかった からです。

理論そのものは頻度主義統計学よりも前にありましたが，「分布を生み出す」には相当な計算量が必要でした。昨今はPCの発展により，一般的なパソコンのパワーでも自由自在に分布を作り出せるようになってきました。

ソフトウェア上でも発展がありました。MCMC法という乱数発生メカニズムを応用すれば，自由に分布を作ることができることがわかり，さらにNUTアルゴリズムをいれると自分の欲しい分布を効率よく得ることができることがわかってきました。このNUTアルゴリズムを実装したstanというソフトがRのパッケージ(rstan)としてCRANに登録され，いつものようにパッケージの導入だけで利用できるようになってきたのは，ごく最近のことです。

rstanを使ってみよう

rstanを使うには，Cコンパイラ というコンピュータ言語が必要です。WindowsのひとはRの他にRtoolsをインストールする必要があります。MacのひとはXcodeが必要です。

これらはそれほど問題ないのですが，それに加えてrstan独特の「モデルの表現方法」を習得する必要があります。一歩一歩，確実に進んでいきましょう。

世界一簡単なrstanコード

stanのコードの例です。lavaanのモデルのように，複数の行に渡っていますので，ダブルクォーテーションでくくります。

stancode1 <- "
data{
  int<lower=0> N;
  real Y[N];
}

parameters{
  real mu;
  real<lower=0> s;
}

model{
  Y ~ normal(mu,s);
  s ~ cauchy(0,5);
}
"

これだけではなにも生じませんね。stancode1オブジェクトにstanのコードを書き込んだだけだからです。データの分析例は次のようになります。

世界一簡単なコードを走らせてみよう

library(rstan)

## Loading required package: ggplot2
## rstan (Version 2.8.2, packaged: 2015-11-26 15:27:02 UTC, GitRev: 05c3d0058b6a)
## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())

# 起動後の呪文
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# 仮データ作成
n <- 100
y <- rnorm(100,50,10)

datastan <- list(Y=y,N=n)
fit <- stan(model_code=stancode1,data=datastan, iter=1000, chain=4)
fit

## Inference for Stan model: 530bc7127c77013eede69849a7cd726e.
## 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
## mu     48.96    0.03 1.03   46.98   48.24   48.95   49.67   50.98  1158
## s      10.49    0.02 0.75    9.15    9.98   10.45   10.97   12.09  1261
## lp__ -284.57    0.04 1.00 -287.44 -284.98 -284.26 -283.87 -283.60   781
##      Rhat
## mu      1
## s       1
## lp__    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Dec 16 09:54:57 2015.
## 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).

うまく推定されているでしょうか。コードの意味を見ながら，いろいろ考えてみてください。

グラフィカルな結果を見たい場合は，shinystanパッケージを使うのがオススメです。

library(shinystan)
launch_shinystan(fit)

世界一簡単なコードをハックしてみよう

最初の例では，平均50，標準偏差10の正規乱数を発生させて，そのパラメタが復元するかどうかを確認しただけでした。実際のデータに使うにはどうしたらいいのでしょう？

例えば，野球選手の身長のデータを見てみましょう。

b.data$height

##  [1] 179 177 178 176 173 173 183 178 174 178 192 185 183 177 171 175 191
## [18] 180 180 176 177 175 179 182 186 180 185 183 171 188 180 174 177 180
## [35] 176 175 177 178 185 177 173 196 175 180 179 188 185 183 175 182 176
## [52] 177 180 180 174 178 182 198 183 180 180 187 182 183 188 180 186 177
## [69] 188 175 185 185 177 187 180 175 194 180 185 180

mean(b.data$height)

## [1] 180.525

sd(b.data$height)

## [1] 5.568844

hist(b.data$height,breaks=20)

普通の人よりは大きいですが，正規分布に従っているようです。このデータから，野球選手の身長はどういう分布をするか，stanで生成してみましょう。

datastan <- list(Y=b.data$height,N=length(b.data$height))
fit.height <- stan( model_code=stancode1,data=datastan, iter=1000, chain=4)
fit.height

## Inference for Stan model: 530bc7127c77013eede69849a7cd726e.
## 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
## mu    180.56    0.02 0.61  179.34  180.17  180.57  180.98  181.74  1352
## s       5.61    0.01 0.43    4.83    5.30    5.58    5.89    6.52  1230
## lp__ -176.92    0.04 0.93 -179.44 -177.31 -176.62 -176.25 -175.99   655
##      Rhat
## mu      1
## s       1
## lp__    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Dec 16 09:55:33 2015.
## 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).

「記述統計量と何も変わらないんじゃ・・・」と思うかもれませんが，生成された分布のパラメータが分布している(点じゃない)ところがベイジアンです。

fit.height.EXT <- extract(fit.height)

# 平均が分布する
fit.height.EXT$mu 
hist(fit.height.EXT$mu,breaks=20)

# 標準偏差が分布する
fit.height.EXT$s
hist(fit.height.EXT$s,breaks=20)

分布と確率で遊ぶ

分布の形を変えてみる

正規分布以外の分布の形を推定してみましょう。例えば，ホームランの本数はポアソン分布に従いますので，ポアソン分布のパラメータを推定してみます。ポアソン分布は，パラメータが一つ(正の実数)だけです。

stancode2 <- "
data{
  int<lower=0> N;
  int Y[N];
}

parameters{
  real<lower=0> lambda;
}

model{
  Y ~ poisson(lambda);
}
"

このコードを走らせてみると，

datastan <- list(Y=b.data$HR,N=length(b.data$HR))
fit.HR <- stan( model_code=stancode2,data=datastan, iter=1000, chain=4)
fit.HR

## Inference for Stan model: 5a5b61809a17f1d07d23efd05b1bef41.
## 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
## lambda   10.46    0.02 0.37    9.76   10.20   10.45   10.71   11.22   502
## lp__   1129.93    0.03 0.75 1127.73 1129.79 1130.21 1130.40 1130.45   700
##        Rhat
## lambda 1.01
## lp__   1.01
## 
## Samples were drawn using NUTS(diag_e) at Wed Dec 16 09:56:09 2015.
## 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).

平均的な本数がわかります。

確率的な問いを立ててみる

さて，ここでメジャーのバッターは年間平均11本ホームランを打つことがわかっているとします。日本の野球選手たちがメジャーの選手に勝つ割合（平均的なHR本数が増える割合)はどれぐらいでしょうか？

このような問いは，頻度主義的には答えにくいのですが，分布を発生させているベイジアンであれば，簡単に答えることができます。コードを次のように，少し書き加えます。

stancode3 <- "
data{
  int<lower=0> N;
  int Y[N];
}

parameters{
  real<lower=0> lambda;
}

model{
  Y ~ poisson(lambda);
}

generated quantities{
  real<lower=0> p;
  p <- step(lambda-11);
}
"

このgenerated quantities（生成量）のところで，平均\(\lambda\)が11よりも大きくなった割合を求めています。このコードを走らせてみましょう。

datastan <- list(Y=b.data$HR,N=length(b.data$HR))
fit.HR2 <- stan( model_code=stancode3,data=datastan, iter=1000, chain=4)
fit.HR2

## Inference for Stan model: c657ae0d897de358539d766549ea0fc5.
## 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
## lambda   10.48    0.01 0.36    9.79   10.23   10.48   10.71   11.21   760
## p         0.08    0.01 0.27    0.00    0.00    0.00    0.00    1.00  1066
## lp__   1129.96    0.02 0.70 1127.85 1129.78 1130.23 1130.41 1130.45   881
##        Rhat
## lambda 1.00
## p      1.00
## lp__   1.01
## 
## Samples were drawn using NUTS(diag_e) at Wed Dec 16 09:56:44 2015.
## 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).

このように，分布の形を変えることができ，また仮説の立て方も平均値や分散，分位点などに対して考えられるので，柔軟なモデリングができます。

検定を考えてみる

例えば，二群の平均値の差の検定はt検定でした。セ・リーグとパ・リーグの選手の平均身長に差があるかどうかを考えるときに，

二群の母平均に差がない，という帰無仮説を立てる
有意水準を考える（普通は5％）
t検定量を算出する

という流れでしたね。分析方法としては，次のように行います。

t.test(height~CP,data=b.data)

## 
##  Welch Two Sample t-test
## 
## data:  height by CP
## t = 0.010512, df = 65.106, p-value = 0.9916
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  -2.519807  2.546474
## sample estimates:
## mean in group C mean in group P 
##        180.5333        180.5200

これをベイジアンでモデリングするとどうすればよいでしょうか。

二つの群が正規分布から発生している
二つの群から作った正規分布の平均の差を・・・みる。

やってみます。

stancode4 <- "
data{
  int<lower=0> N1;
  int<lower=0> N2;
  real Y1[N1];
  real Y2[N2];
}

parameters{
  real mu1;
  real<lower=0> s1;
  real mu2;
  real<lower=0> s2;
}

model{
  Y1 ~ normal(mu1,s1);
  Y2 ~ normal(mu2,s2);
  s1 ~ cauchy(0,5);
  s2 ~ cauchy(0,5);
}

generated quantities{
  real diff;
  diff <- mu1- mu2;
}
"

C <- subset(b.data,b.data$CP=="C")
P <- subset(b.data,b.data$CP=="P")

datastan <- list(N1=nrow(C),N2=nrow(P),Y1=C$height,Y2=P$height)
fit.height3 <- stan(model_code=stancode4,data=datastan, iter=1000, chain=4)
fit.height3

## Inference for Stan model: f99622f4a01f1f6007214ef3069c5d73.
## 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
## mu1   180.54    0.03 1.00  178.53  179.90  180.53  181.19  182.50  1192
## s1      5.48    0.02 0.74    4.25    4.96    5.39    5.93    7.12  1005
## mu2   180.50    0.02 0.80  178.93  179.95  180.50  181.00  182.11  1393
## s2      5.85    0.02 0.57    4.84    5.45    5.79    6.20    7.07  1213
## diff    0.05    0.04 1.29   -2.52   -0.81    0.07    0.90    2.57  1178
## lp__ -176.88    0.05 1.37 -180.24 -177.65 -176.58 -175.82 -175.12   696
##      Rhat
## mu1     1
## s1      1
## mu2     1
## s2      1
## diff    1
## lp__    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Dec 16 09:57:19 2015.
## 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).

「差がないという帰無仮説が棄却できない」といった言い回しではなく，差があるとしても平均これぐらい（だから考えるほどじゃない）というような言い方ができるかと思います。この差の出し方を工夫することで，もっと色々な仮説を考えることができます。例えば，

セリーグの平均身長がパリーグの平均身長を上回る確率はどれぐらいあるのか
セリーグとパリーグの平均身長が1cm以上離れる確率はどれぐらいあるのか

といった話から，分布の形を考慮しつつ，

セリーグの平均年収がパリーグの平均年収を上回る確率はどれぐらいあるのか
セリーグのHR本数の平均がパリーグのHR本数の平均を3本上回る確率は

といったことを検証することができます。

例；セリーグの平均年収がパリーグの平均年収を上回る確率はどれぐらいあるのか

stancode5 <- "
data{
  int<lower=0> N1;
  int<lower=0> N2;
  real Pay1[N1];
  real Pay2[N2];
}

parameters{
  real<lower=0> mu1;
  real<lower=0> s1;
  real<lower=0> mu2;
  real<lower=0> s2;
}

model{
  Pay1 ~ lognormal(mu1,s1);
  Pay2 ~ lognormal(mu2,s2);
  s1 ~ cauchy(0,5);
  s2 ~ cauchy(0,5);
}

generated quantities{
  real delta;
  real prob;
  delta <- mu1- mu2;
  prob <- step(delta);
}
"

datastan <- list(N1=nrow(C),N2=nrow(P),Pay1=C$pay,Pay2=P$pay)
fit.pay <- stan(model_code=stancode5,data=datastan, iter=1000, chain=4)
fit.pay

## Inference for Stan model: fbd201fdb11a1f7c685200291379a943.
## 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
## mu1    2.53    0.00 0.12  2.29  2.45  2.53  2.61  2.78  1220 1.00
## s1     0.63    0.00 0.09  0.50  0.57  0.62  0.68  0.84  1024 1.00
## mu2    2.49    0.00 0.09  2.31  2.43  2.49  2.55  2.66  1293 1.00
## s2     0.64    0.00 0.07  0.52  0.59  0.63  0.68  0.79   915 1.00
## delta  0.05    0.00 0.15 -0.24 -0.05  0.05  0.15  0.33  1322 1.00
## prob   0.62    0.01 0.49  0.00  0.00  1.00  1.00  1.00  1500 1.00
## lp__  -1.46    0.06 1.52 -5.31 -2.19 -1.08 -0.38  0.37   651 1.01
## 
## Samples were drawn using NUTS(diag_e) at Wed Dec 16 09:57:55 2015.
## 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).

本日の課題

バッターのデータを使い，平均盗塁数を推定してください。盗塁データはstealです。
バッターのデータを使い，セリーグの平均盗塁数がパリーグの平均盗塁数を上回る確率を求めてください

Lesson9.Bayesian Modeling

Kosugitti

2015年12月18日