Lesson8.Next Step for Bayesian Modeloing
Kosugitti
2017年12月15日
毎回やる儀式
- インターネットに接続しておいてください（このサイトが見えていればOK）。
- Rstudioを起動し，この授業のプロジェクトを開いてください。(File>Open Project)
- File>New File>R scriptで新しいコードファイルを作ってください。
謝罪
先週は体調不良で休講しちゃってごめんね。
本日のお題
昨年の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パーセンタイルを参考に答えてください。
- 評価のばらつきが最も大きい評定者は誰ですか。中央値を参考にして（行番号で）答えてください。
- 評価のばらつきが最も小さい評定者は誰ですか。中央値を参考にして（行番号で）答えてください。
data{
int<lower=1> N; //number of players
int<lower=1> M; //number of rators
real X[M,N]; // scores
}
parameters{
real<lower=0,upper=100> theta[N];
real<lower=0> sig[M];
}
model{
for(n in 1:N){
for(m in 1:M){
X[m,n] ~ normal(theta[n],sig[m]);
}
}
theta ~ normal(50,10);
sig ~ cauchy(0,5);
}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())dat <- read.csv("owarai1.csv",fileEncoding = "UTF-8")[,-1]
datastan <- list(N=ncol(dat),M=nrow(dat),X=dat)
fit <- sampling(owarai.stan,datastan)##
## SAMPLING FOR MODEL '281c2d42b5a8d5ec0880bb9f57bf41de' NOW (CHAIN 1).
##
## Gradient evaluation took 2.6e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.26 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 0.157307 seconds (Warm-up)
## 0.138783 seconds (Sampling)
## 0.29609 seconds (Total)
##
##
## SAMPLING FOR MODEL '281c2d42b5a8d5ec0880bb9f57bf41de' NOW (CHAIN 2).
##
## Gradient evaluation took 1.4e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 0.156523 seconds (Warm-up)
## 0.183265 seconds (Sampling)
## 0.339788 seconds (Total)
##
##
## SAMPLING FOR MODEL '281c2d42b5a8d5ec0880bb9f57bf41de' NOW (CHAIN 3).
##
## Gradient evaluation took 1.3e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.13 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 0.158337 seconds (Warm-up)
## 0.124903 seconds (Sampling)
## 0.28324 seconds (Total)
##
##
## SAMPLING FOR MODEL '281c2d42b5a8d5ec0880bb9f57bf41de' NOW (CHAIN 4).
##
## Gradient evaluation took 1.4e-05 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 0.14 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 200 / 2000 [ 10%] (Warmup)
## Iteration: 400 / 2000 [ 20%] (Warmup)
## Iteration: 600 / 2000 [ 30%] (Warmup)
## Iteration: 800 / 2000 [ 40%] (Warmup)
## Iteration: 1000 / 2000 [ 50%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 1200 / 2000 [ 60%] (Sampling)
## Iteration: 1400 / 2000 [ 70%] (Sampling)
## Iteration: 1600 / 2000 [ 80%] (Sampling)
## Iteration: 1800 / 2000 [ 90%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 0.157019 seconds (Warm-up)
## 0.188518 seconds (Sampling)
## 0.345537 seconds (Total)print(fit)## Inference for Stan model: 281c2d42b5a8d5ec0880bb9f57bf41de.
## 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%
## theta[1] 81.86 0.06 2.34 77.34 80.23 81.85 83.50 86.33
## theta[2] 77.06 0.07 2.54 71.80 75.37 77.23 78.84 81.57
## theta[3] 81.62 0.06 2.28 77.00 80.12 81.67 83.16 86.00
## theta[4] 74.64 0.05 2.27 70.12 73.11 74.73 76.22 78.94
## theta[5] 83.32 0.04 2.00 79.14 82.02 83.39 84.69 87.04
## theta[6] 83.95 0.05 2.18 79.65 82.49 83.94 85.39 88.28
## theta[7] 80.57 0.05 2.24 75.84 79.16 80.77 82.11 84.58
## theta[8] 85.74 0.04 2.04 81.60 84.41 85.81 87.16 89.52
## theta[9] 83.15 0.03 1.94 79.39 81.82 83.11 84.47 87.00
## sig[1] 6.17 0.04 1.90 3.42 4.79 5.84 7.19 10.58
## sig[2] 13.39 0.07 3.57 8.33 10.93 12.77 15.09 22.13
## sig[3] 19.97 0.09 4.91 12.73 16.45 19.23 22.58 31.45
## sig[4] 10.30 0.07 2.99 6.11 8.23 9.78 11.73 17.59
## sig[5] 8.36 0.06 2.43 4.93 6.67 7.95 9.57 14.12
## sig[6] 4.01 0.05 1.55 1.56 2.92 3.82 4.87 7.54
## sig[7] 7.96 0.04 2.12 4.95 6.47 7.59 9.07 13.11
## sig[8] 9.04 0.05 2.53 5.51 7.30 8.60 10.27 15.00
## sig[9] 4.65 0.05 1.66 2.00 3.54 4.44 5.54 8.45
## sig[10] 7.31 0.05 2.11 4.21 5.86 6.99 8.37 12.30
## sig[11] 9.15 0.05 2.57 5.49 7.37 8.70 10.44 15.43
## sig[12] 7.05 0.04 1.89 4.35 5.73 6.72 8.02 11.86
## lp__ -290.94 0.08 3.35 -298.42 -293.02 -290.61 -288.52 -285.31
## n_eff Rhat
## theta[1] 1442 1
## theta[2] 1359 1
## theta[3] 1534 1
## theta[4] 1715 1
## theta[5] 3118 1
## theta[6] 1804 1
## theta[7] 1718 1
## theta[8] 2600 1
## theta[9] 3502 1
## sig[1] 1931 1
## sig[2] 2463 1
## sig[3] 3169 1
## sig[4] 1892 1
## sig[5] 1816 1
## sig[6] 1136 1
## sig[7] 3324 1
## sig[8] 2629 1
## sig[9] 1281 1
## sig[10] 1973 1
## sig[11] 2372 1
## sig[12] 2690 1
## lp__ 1686 1
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 10:23:16 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).# 平均的な実力が最も高い芸人は和牛
# 最も高い実力をもっている可能性がある芸人も和牛
# 低い可能性はスリムクラブ
# 大きいのは久保山，小さいのは浜崎自由なモデリングの世界
以下三つの記事を読んでください。
参考記事1;最強のM-1漫才師は誰だ 参考記事2;「最強のM-1漫才師は誰だ」へのチャレンジ 参考記事3;「最強のM-1漫才師は誰だ」シリーズへの挑戦
自分ならこういうことを考えたい，というのを語り合いましょう
検定とモデリング(講義)
心理学の研究では，t検定や分散分析といった，「平均値を比べる」ということがよく行われます。 また，回帰分析によって，独立変数の従属変数に対する影響力を検証したりします。 この二つは違う研究法なのでしょうか？ いいえ，いずれもモデリングという意味では同じなのです。特に，平均値を比べるモデリングのことを一般線形モデルといいます。
平均値の差を検証するモデリング
今年のM-1グランプリの評定値をつかって，次の「平均値」を比べたいと思います。
- 和牛ととろサーモン
- ミキと和牛ととろサーモン
m1 <- read.csv("M1score.csv",fileEncoding = "UTF-8",na.strings = ".")
m1.13 <- subset(m1,m1$年代==13)
m1.miki <- as.numeric(m1.13[3,3:9]) # データの3行目，3列目ー9列目の数字
m1.toro <- as.numeric(m1.13[8,3:9]) # データの8行目，3列目ー9列目の数字
m1.wagu <- as.numeric(m1.13[2,3:9]) # データの2行目，3列目ー9列目の数字2群を比べるモデリング
次のコードをファイル名「2groups.stan」として保存してください。
data{
int N1;
int N2;
real X1[N1];
real X2[N2];
}
parameters{
real<lower=0,upper=100> mu;
real diff;
real<lower=0> sigma1;
real<lower=0> sigma2;
}
model{
//LIKELIHOOD
X1 ~ normal(mu,sigma1);
X2 ~ normal(mu+diff,sigma2);
//PRIOR
mu ~ uniform(0,100);
sigma1 ~ cauchy(0,5);
sigma2 ~ cauchy(0,5);
diff ~ normal(0,100);
}その上で，次のコードで読み込んで実行してください。
library(rstan)
options(mc.cores = parallel::detectCores())
datastan <- list(N1=7,N2=7,X1=m1.toro,X2=m1.wagu)
model <- stan_model("2groups.stan")
fit <- sampling(model,datastan)
fit## Inference for Stan model: 2groups.
## 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 92.15 0.02 0.92 90.25 91.63 92.18 92.71 93.92 1574 1
## diff 1.11 0.03 1.09 -1.01 0.44 1.10 1.78 3.34 1576 1
## sigma1 2.28 0.02 0.84 1.25 1.72 2.09 2.63 4.43 1549 1
## sigma2 1.41 0.01 0.57 0.75 1.04 1.29 1.64 2.71 1527 1
## lp__ -10.94 0.05 1.66 -14.86 -11.81 -10.60 -9.70 -8.86 1175 1
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 10:23:48 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,pars="diff",show_density=T)## ci_level: 0.8 (80% intervals)## outer_level: 0.95 (95% intervals)
結果について，次のことを考えてみよう。
- とろサーモンと和牛はどっちの評価の方が上か？
- とろサーモンと和牛はどっちの評価の方がブレている？
- とろサーモンと和牛の実力差はどれぐらいか？
- とろサーモンと和牛の実力差は100%確実なものか？
3群を比べるモデリング
次のコードをファイル名「3groups.stan」として保存してください。
data{
int N1;
int N2;
int N3;
real X1[N1];
real X2[N2];
real X3[N3];
}
parameters{
real<lower=0,upper=100> mu;
real d1;
real d2;
real<lower=0> sigma1;
real<lower=0> sigma2;
real<lower=0> sigma3;
}
transformed parameters{
real diff[3];
diff[1] = d1;
diff[2] = d2;
diff[3] = 0 - (d1+d2);
}
model{
//LIKELIHOOD
X1 ~ normal(mu+diff[1],sigma1);
X2 ~ normal(mu+diff[2],sigma2);
X3 ~ normal(mu+diff[3],sigma3);
//PRIOR
mu ~ uniform(0,100);
sigma1 ~ cauchy(0,5);
sigma2 ~ cauchy(0,5);
sigma3 ~ cauchy(0,5);
d1 ~ normal(0,100);
d2 ~ normal(0,100);
}
その上で，次のコードで読み込んで実行してください。
datastan <- list(N1=7,N2=7,N3=7,X1=m1.toro,X2=m1.wagu,X3=m1.miki)
model2 <- stan_model("3groups.stan")
fit2 <- sampling(model2,datastan)
fit2## Inference for Stan model: 3groups.
## 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 92.74 0.01 0.44 91.84 92.48 92.75 93.01 93.63 2328 1
## d1 -0.66 0.01 0.70 -2.06 -1.08 -0.65 -0.22 0.70 2197 1
## d2 0.55 0.01 0.55 -0.52 0.20 0.54 0.91 1.67 2143 1
## sigma1 2.29 0.02 0.83 1.26 1.74 2.09 2.63 4.52 1883 1
## sigma2 1.39 0.01 0.53 0.75 1.04 1.27 1.61 2.70 2461 1
## sigma3 1.94 0.01 0.70 1.07 1.47 1.78 2.22 3.71 2242 1
## diff[1] -0.66 0.01 0.70 -2.06 -1.08 -0.65 -0.22 0.70 2197 1
## diff[2] 0.55 0.01 0.55 -0.52 0.20 0.54 0.91 1.67 2143 1
## diff[3] 0.10 0.01 0.63 -1.13 -0.29 0.10 0.50 1.38 4000 1
## lp__ -17.90 0.06 1.99 -22.69 -18.96 -17.47 -16.45 -15.23 1309 1
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 10:24:21 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).結果について，次のことを考えてみよう。
- とろサーモンと和牛はどっちの評価の方が上か？
- とろサーモンとミキはどっちの評価の方が上か？
- 和牛とミキはどっちの評価の方が上か？
- 審査員の評価ブレが大きかったのはどのコンビ？
- とろサーモンとミキの実力差はどれぐらいか？
- ミキと和牛の実力差は100%確実なものか？
コードを洗練させる
ちなみに，上のコードは次のように書き直すこともできます。データのサイズがわかっているので，読み込む時に\(X[3,7]\)に固定し，同じことを何度も書かなくていいようにします。
data{
real X[3,7];
}
parameters{
real<lower=0,upper=100> mu;
real d[2];
real<lower=0> sigma[3];
}
transformed parameters{
real diff[3];
diff[1] = d[1];
diff[2] = d[2];
diff[3] = 0 - sum(d);
}
model{
//LIKELIHOOD
for(i in 1:3){
X[i] ~ normal(mu+diff[i],sigma[i]);
}
//PRIOR
mu ~ uniform(0,100);
sigma ~ cauchy(0,5);
d ~ normal(0,100);
}データの形が変わったので，RからStanにデータを渡す時の渡し方も変えます。
datX <- m1.13[c(3,8,2),3:9]
datastan <- list(X=datX)
model2b <- stan_model("3groups_b.stan")
fit2b <- sampling(model2b,datastan)結果は変わってないです。モデルとしては同じだから。
print(fit2b)## Inference for Stan model: 3groups_b.
## 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 92.75 0.01 0.45 91.83 92.47 92.77 93.04 93.60 1476 1
## d[1] 0.13 0.02 0.64 -1.16 -0.26 0.12 0.51 1.45 1635 1
## d[2] -0.66 0.02 0.69 -2.07 -1.07 -0.64 -0.22 0.64 1076 1
## sigma[1] 1.96 0.02 0.73 1.06 1.47 1.80 2.27 3.81 1902 1
## sigma[2] 2.29 0.03 0.86 1.29 1.73 2.11 2.64 4.44 996 1
## sigma[3] 1.40 0.01 0.52 0.76 1.05 1.30 1.62 2.74 2350 1
## diff[1] 0.13 0.02 0.64 -1.16 -0.26 0.12 0.51 1.45 1635 1
## diff[2] -0.66 0.02 0.69 -2.07 -1.07 -0.64 -0.22 0.64 1076 1
## diff[3] 0.54 0.01 0.56 -0.56 0.17 0.54 0.89 1.67 2307 1
## lp__ -17.92 0.06 1.99 -22.67 -19.00 -17.56 -16.44 -15.20 1080 1
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 10:24:52 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).データは行列型でなく，縦長に変形する方が便利なことがあります。欠損値が含まれている場合は，縦長にして欠損値をカットするとよいでしょう。その場合，データの何列目が，誰の，どのデータか，というのを指し示す指標(インデックス)変数が必要になります。
縦長にする場合は，tidyrパッケージを使うと便利です。
# 行列型（横長，wide型）
m1.13[c(3,8,2),2:9]## 演者 上沼恵美子 松本人志 博多大吉 春風亭小朝 中川礼二 渡辺正行
## 3 ミキ 95 94 91 92 91 94
## 8 とろサーモン 93 92 93 93 93 93
## 2 和牛 95 93 94 94 93 92
## オール巨人
## 3 93
## 8 88
## 2 92# 縦長型(long型)
library(tidyr)##
## Attaching package: 'tidyr'## The following object is masked from 'package:rstan':
##
## extractm1.13.long <- gather(m1.13[c(3,8,2),2:9],rator,score,-演者)
head(m1.13.long)## 演者 rator score
## 1 ミキ 上沼恵美子 95
## 2 とろサーモン 上沼恵美子 93
## 3 和牛 上沼恵美子 95
## 4 ミキ 松本人志 94
## 5 とろサーモン 松本人志 92
## 6 和牛 松本人志 93# 欠損値を消しておく
m1.13.long <- na.omit(m1.13.long)
# 中身の確認
m1.13.long## 演者 rator score
## 1 ミキ 上沼恵美子 95
## 2 とろサーモン 上沼恵美子 93
## 3 和牛 上沼恵美子 95
## 4 ミキ 松本人志 94
## 5 とろサーモン 松本人志 92
## 6 和牛 松本人志 93
## 7 ミキ 博多大吉 91
## 8 とろサーモン 博多大吉 93
## 9 和牛 博多大吉 94
## 10 ミキ 春風亭小朝 92
## 11 とろサーモン 春風亭小朝 93
## 12 和牛 春風亭小朝 94
## 13 ミキ 中川礼二 91
## 14 とろサーモン 中川礼二 93
## 15 和牛 中川礼二 93
## 16 ミキ 渡辺正行 94
## 17 とろサーモン 渡辺正行 93
## 18 和牛 渡辺正行 92
## 19 ミキ オール巨人 93
## 20 とろサーモン オール巨人 88
## 21 和牛 オール巨人 92gather関数の挙動をよく理解しておきましょう。基本的には全部縦長にする，されたくない変数はマイナスの記号でどけておく，というものです。
そしてこれを使ってモデリングする場合は，「演者」や「rator」を指標変数＝数字にして渡してやる必要があります。今は文字として入っていますから，factor型にして数字とラベルの対応をつけ，そこの数字だけ抜き出す，という二段階を経ます。
# factor型にする
factor(m1.13.long$演者)## [1] ミキ とろサーモン 和牛 ミキ とろサーモン
## [6] 和牛 ミキ とろサーモン 和牛 ミキ
## [11] とろサーモン 和牛 ミキ とろサーモン 和牛
## [16] ミキ とろサーモン 和牛 ミキ とろサーモン
## [21] 和牛
## Levels: とろサーモン ミキ 和牛# factor型の数字だけ抜き出す
as.numeric(factor(m1.13.long$演者))## [1] 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3 2 1 3これらの指標変数を使ったモデリングは次のようになります。ファイル名を3group_c.stanにして保存してください。
data{
int<lower=0> L; //データの長さ
int playerID[L]; // l 行目のプレイヤー・インデックス
real X[L]; // l行目のスコア
}
parameters{
real<lower=0,upper=100> mu;
real d[2];
real<lower=0> sigma[3];
}
transformed parameters{
real diff[3];
diff[1] = d[1];
diff[2] = d[2];
diff[3] = 0 - sum(d);
}
model{
//LIKELIHOOD
for(l in 1:L){
X[l] ~ normal(mu+diff[playerID[l]],sigma[playerID[l]]);
}
//PRIOR
mu ~ uniform(0,100);
sigma ~ cauchy(0,5);
d ~ normal(0,100);
}データの渡し方に注意し，実行して見ます。結果が変わらないことにも注目です。
datastan <- list(L=nrow(m1.13.long),playerID=as.numeric(factor(m1.13.long$演者)),X=m1.13.long$score)
model2c <- stan_model("3groups_c.stan")
fit2c <- sampling(model2c,datastan)
fit2c## Inference for Stan model: 3groups_c.
## 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 92.76 0.01 0.43 91.88 92.50 92.75 93.04 93.62 2631 1
## d[1] -0.63 0.02 0.68 -2.04 -1.04 -0.62 -0.20 0.72 1793 1
## d[2] 0.10 0.01 0.62 -1.11 -0.29 0.09 0.48 1.30 1978 1
## sigma[1] 2.26 0.02 0.78 1.27 1.73 2.08 2.60 4.29 1883 1
## sigma[2] 1.94 0.02 0.74 1.07 1.47 1.78 2.23 3.70 1791 1
## sigma[3] 1.42 0.01 0.53 0.76 1.06 1.31 1.63 2.81 1741 1
## diff[1] -0.63 0.02 0.68 -2.04 -1.04 -0.62 -0.20 0.72 1793 1
## diff[2] 0.10 0.01 0.62 -1.11 -0.29 0.09 0.48 1.30 1978 1
## diff[3] 0.53 0.01 0.54 -0.51 0.18 0.53 0.86 1.65 4000 1
## lp__ -17.86 0.05 1.94 -22.68 -18.92 -17.48 -16.45 -15.19 1297 1
##
## Samples were drawn using NUTS(diag_e) at Thu Dec 14 10:25:25 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にでた10組，全ての実力差を検証してください。一番実力差があったのは，どことどこのコンビのコンビネーションでしょう？ コンビ\(j\)のスコア\(X_j\)が，真の実力\(\mu_j\)をもっており，また評定者\(k\)による分散\(\sigma_k\)も考慮して，\(X_{jk} \sim N(\mu_j,\sigma_k)\)と考えます。 プレイヤーは10組で，各組みの平均\(\mu_j\)は\(\mu_j = \mu + diff_j\)と考えてモデリングしましょう。この時，\(\sum_{j=1}^{10} diff_j = 0\)とする(全体平均からの偏差とする)ことがポイントです。