毎回やる儀式

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

今回必要なパッケージ

前期の授業で導入しているかと思いますが，アップデートがあるかもしれませんので要注意。次のコードを走らせて，分析結果が出るようならオーケーです。試して見てください。

library(rstan)

## Loading required package: ggplot2

## Loading required package: StanHeaders

## rstan (Version 2.16.2, packaged: 2017-07-03 09:24:58 UTC, GitRev: 2e1f913d3ca3)

## 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())

stancode <- '
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 {
  target += normal_lpdf(eta | 0, 1);
  target += normal_lpdf(y | theta, sigma);
}
'

schools_dat <- list(J = 8, 
                    y = c(28,  8, -3,  7, -1,  1, 18, 12),
                    sigma = c(15, 10, 16, 11,  9, 11, 10, 18))

model.ex <- stan_model(model_code=stancode,model_name="school")
fit.samp <- sampling(model.ex, data = schools_dat, iter = 1000, chains = 2)
fit.samp

## Inference for Stan model: school.
## 2 chains, each with iter=1000; warmup=500; thin=1; 
## post-warmup draws per chain=500, total post-warmup draws=1000.
## 
##            mean se_mean   sd   2.5%    25%    50%    75%  97.5% n_eff Rhat
## mu         7.89    0.26 5.08  -2.43   4.57   7.87  11.31  18.25   394 1.01
## tau        6.35    0.25 5.22   0.28   2.48   5.04   9.03  19.81   425 1.00
## eta[1]     0.39    0.03 0.95  -1.60  -0.22   0.43   0.98   2.20  1000 1.00
## eta[2]    -0.08    0.03 0.89  -1.80  -0.67  -0.10   0.49   1.77  1000 1.00
## eta[3]    -0.18    0.03 0.88  -1.93  -0.75  -0.18   0.43   1.48  1000 1.00
## eta[4]     0.03    0.03 0.93  -1.78  -0.57   0.00   0.60   1.95  1000 1.00
## eta[5]    -0.39    0.03 0.87  -2.10  -0.97  -0.37   0.19   1.34  1000 1.00
## eta[6]    -0.23    0.03 0.91  -2.13  -0.79  -0.25   0.36   1.60   904 1.00
## eta[7]     0.31    0.03 0.90  -1.48  -0.30   0.31   0.91   2.04  1000 1.00
## eta[8]     0.05    0.03 0.96  -1.78  -0.62   0.05   0.68   1.95  1000 1.00
## theta[1]  11.26    0.25 8.00  -2.25   6.02  10.40  15.63  30.41  1000 1.00
## theta[2]   7.59    0.19 6.14  -4.81   3.74   7.59  11.63  19.58  1000 1.00
## theta[3]   6.19    0.25 7.78 -10.92   1.94   6.86  11.07  20.47  1000 1.00
## theta[4]   7.62    0.22 6.87  -5.94   3.33   7.91  12.25  20.35  1000 1.00
## theta[5]   5.06    0.21 6.74 -10.32   1.55   5.64   9.32  17.01  1000 1.00
## theta[6]   6.06    0.22 6.85  -9.83   2.07   6.54  10.42  19.15  1000 1.00
## theta[7]  10.54    0.22 6.86  -1.75   5.84  10.02  14.55  25.39  1000 1.00
## theta[8]   8.33    0.25 7.88  -7.97   3.64   8.02  12.72  25.88  1000 1.00
## lp__     -39.68    0.15 2.60 -45.59 -41.21 -39.42 -37.84 -35.31   304 1.01
## 
## Samples were drawn using NUTS(diag_e) at Wed Nov 29 09:59:51 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).

ベイジアンになるための二つの入り口

ベイズ統計学を使ってデータを分析するのには，二つの入り口があります。

従来通りの推定方法ではうまく行かないので，推定法としてのベイズ

この授業の性質的にはこちらが王道です。

今までもデータに対してモデルを当てはめる，ということはやって来ました。回帰分析も，因子分析も，構造方程式モデリングも，データに対してモデルを当てはめる，という考え方です。だから適合度指標(調整済み\(R^2\)，因子数の決定基準(BFIとか)，適合度指標(CFI))という考え方があるわけで。

当てはめる方法は色々あって，最も単純なものは「最小二乗法」というやりかた。誤差を一番少なくするという考え方。二つ目は，「最尤法」というやりかた。これは統計学的な理由が背後にあります。

この二つの方法があれば，大概のモデルはうまく行ってたのですが，もっと色々やりたいとなると「この方法じゃうまく推定できないよ」「推定したけど理論的にありえない数字が出ちゃったよ」というエラー，警告がでたのでした。これを解決するために別の技術が必要，というところで出て来たのがベイズ推定法，というものです。

ちなみに，ベイズ推定ができるパッケージも色々開発されていますから，いちいちゼロからコードを書かなくてもベイズを使うことはできるんですよ。

データ生成モデルを考えるルート

前期の授業でベイズの話は触れて来ましたが，そのとき考えたのは「このデータを作るメカニズムはどうなっているか」ということでした。このルートから考えて，メカニズムの形を探って行くというのが「モデリング」の発想です。皆さんはこちらの考え方にも慣れておいて欲しいと思います。

ルートが違えばマナーも違う

前者の考え方でいくと，「いままでの方法がダメだから新しい方法を使うよ」というだけの話なのですが，この新しい方法を使うには考え方の方も変えなければならないのです。従来の分析法を使っている人にとっては，これがハードルになります（が，皆さんは初めからベイジアンでいくので，それほど抵抗を感じることはないと思います）。

考えなければならない点とは，

点推定ではなく分布で考えなければならない
事前分布というのを置かなければならない
事後分布を得るのはすごく難しいので，乱数発生法の力を借りて近似することになる

というところでしょうか。全体を通じて分布で考えるというのは，「真実は常に一つ！」という発想ではなく「この辺に答えがあるっぽい！」なので，従来型の考え方に慣れ親しんで来た人にはちょっと頭を切り替える必要があります。そもそも，ベイジアンが答えようとしているのは「データから何を，どの程度信用するべきか」という話であって，「データから何をなすべきかを示す」というような考え方とはちょっと違うのです。

問いの形が違うから答え方も違う，と思っていただいても結構。

とはいえ，本講義ではそこまで悩まなくても結構です。とりあえず何ができるかやってみましょう。

怖い本を例に

RQ.7人の科学者

実験スキルが大きく異なる7人の科学者が，全員同じ量について測定を行う。彼らが得た数値結果は次の通り。

x <- c(-27.020,3.570,8.191,9.808,9.603,9.945,10.056)

直感的には，最初の二人の科学者はひどく適性を欠いた測定者であり，この量の真の値はおそらく10をわずかに下回るくらいであるように思えるのだが・・・？？

モデルのイメージはこんな感じです。モデルイメージRQ1

モデル

data{
  real X[7];
}

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

model{
  for(i in 1:7){
    X[i] ~ normal(mu,sigma[i]);
    sigma[i] ~ inv_gamma(0.001,0.001);
  }
  mu ~ normal(0,1000);
}

上のコードを実行するには，ファイル名をつけて保存する必要があります。その際，ファイル名の最後に「.stan」をつけておくのを忘れないように。そうすることで，強調表示がなされミススペルが発見しやすくなります。

上のコードはStanファイルを別に作って分析した方がやりやすいです。これをRの方からRStanパッケージで呼び出します。

プロセス

次のコードではsampling関数で事後分布からのサンプリングを得ていますが，model1.stanのところは各自の命名に合わせて，またファイル名はダブルクォーテーションで囲って実行してください。

fit <- sampling(model1.stan,data=list(X=x))

## Warning: There were 690 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
## http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmup

## Warning: Examine the pairs() plot to diagnose sampling problems

print(fit)

## Inference for Stan model: 923cd487e470d3e29f8dcc09cfbb83e0.
## 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
## mu         9.86    0.01    0.18   9.58  9.80  9.90   9.95   10.07   181
## sigma[1] 222.59   32.31 1667.79  17.10 35.60 68.62 172.96 1070.51  2664
## sigma[2]  46.83   16.56  309.26   2.81  5.17  7.82  16.56  174.28   349
## sigma[3]   6.64    0.79   25.86   0.69  1.29  2.01   4.46   41.04  1063
## sigma[4]   1.13    0.31   11.36   0.00  0.07  0.19   0.39    4.05  1357
## sigma[5]   1.43    0.30    9.14   0.01  0.18  0.40   0.86    6.43   900
## sigma[6]   0.90    0.23    6.40   0.00  0.03  0.12   0.33    4.55   788
## sigma[7]   0.91    0.12    4.33   0.01  0.13  0.24   0.55    5.64  1360
## lp__      -4.60    0.38    2.82 -10.89 -6.31 -4.35  -2.48   -0.21    56
##          Rhat
## mu       1.02
## sigma[1] 1.00
## sigma[2] 1.01
## sigma[3] 1.01
## sigma[4] 1.00
## sigma[5] 1.01
## sigma[6] 1.00
## sigma[7] 1.00
## lp__     1.06
## 
## Samples were drawn using NUTS(diag_e) at Wed Nov 29 10:00: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).

plot(fit,show_density=T,pars=c("mu","sigma"),ci_level=0.5)

## ci_level: 0.5 (50% intervals)

## outer_level: 0.95 (95% intervals)

実行結果から何が読み取れるか，EAP,MEDなどの推定値やHDIなどで表現できることを復習しておきましょう。

RQ. IQの繰り返し測定

ある人々のIQを推定する方法を考える。それぞれの人が複数回のIQテストを受けている。データは\(i=1,2,\cdots, n\)人の人と，\(j=1,2,\cdots,m\)回のくり返しテスト得点についての測度\(X_{ij}\)である。くり返し測定されるテスト得点間の差は，平均0で標準偏差は未知の正規分布と考えられる（誤差だから）。この標準偏差は人を通じて同じだろう（テストの特性だから）。それぞれの人のIQは違っているだろう。

この考えをモデルにします。

モデルのイメージはこんな感じです。モデルイメージRQ2

コードで表現すると次のようになります。

data{
  int N; // number of subjects
  int M; // number of scores
  real X[N,M]; // data
}
parameters{
  real<lower=0> mu[N];
  real<lower=0> sigma;
}
model{
  //model
  for(i in 1:N){
    for(j in 1:M){
      X[i,j] ~ normal(mu[i],sigma);
    }
  }
  //prior
  sigma ~ inv_gamma(0.001,0.001);
  for(i in 1:N){
    mu[i] ~ uniform(0,300);
  }
}

実行してみましょう。

# IQ test
X <- matrix(c(90,95,100,
              105,110,115,
              150,155,160),ncol=3,byrow=T)
N <- 3
M <- 3
fit2 <- sampling(model2.stan,data=list(N=N,M=M,X=X))
print(fit2)

## Inference for Stan model: ba65cf042880382e1fc7126527e862d6.
## 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[1]  95.01    0.06 3.48  88.13  92.99  95.06  97.02 101.92  2948    1
## mu[2] 110.04    0.07 3.51 103.49 107.98 109.99 112.02 117.10  2775    1
## mu[3] 155.01    0.07 3.45 148.31 153.09 154.95 156.93 162.01  2712    1
## sigma   5.75    0.06 2.11   3.20   4.33   5.29   6.60  11.15  1431    1
## lp__   -5.42    0.06 1.88 -10.32  -6.28  -4.97  -4.09  -3.14  1019    1
## 
## Samples were drawn using NUTS(diag_e) at Wed Nov 29 10:00:58 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).

結果をどのように読み取るか，数値の結果を言葉で表現できるように練習しましょう。

本日の課題

昨年のM-1について，皆さんに面白さ評定をしてもらいました。各評定値は，お笑い芸人の面白さの実力に加えて，評定者の癖，特徴も反映していると考えられます。

ここで，お笑い芸人の実力を\(\theta_i\)とします。評定者の特徴\(\sigma_j\)は，評定者ごとにスコアがブレるブレの大きさ，つまりRQ1のような評定の誤差を表していると考えられます。

芸人\(i\)に対する評定者\(j\)の評定値を\(X_{ij}\)とすると，\(X_{ij} \sim N(\theta_i,\sigma_j)\)というモデルが考えられます。このモデルを，owarai1.csvファイルのデータで検証し，次の問いに答えてください。

ただし，事前分布として，\(\theta_i \sim N(50,10)\),\(\sigma_j \sim cauchy(0,5)\)とします。

平均的な実力が最も高い芸人は誰ですか
最も高い実力をもっている可能性がある芸人はだれですか。実力の97.5パーセンタイルを参考に答えてください。
最も点数が低くなる可能性がある芸人は誰ですか。実力の2.5パーセンタイルを参考に答えてください。
評価のばらつきが最も大きい評定者は誰ですか。中央値を参考にして（行番号で）答えてください。
評価のばらつきが最も小さい評定者は誰ですか。中央値を参考にして（行番号で）答えてください。

Lesson6.Introduction to Bayesian Modeloing

Kosugitti

2017年12月01日