![モデル
※[“”]はちゃんとコピペできないので書き直してね
stancode<- “
data {
int<lower=1> n;
int<lower=0> k;
}
parameters {
real<lower=0,upper=1> theta;
}
model {
theta ~ beta(1, 1);
// Observed Counts
k ~ binomial(n, theta);
}”](https://image.slidesharecdn.com/rbayse2015710-150727085101-lva1-app6891/85/R-1-BCM-22-320.jpg)
![Code：長いから覚悟しろよ！
※ “” はコピペすると以下略
? model <- " // Pearson Correlation
data {
int<lower=0> n;
vector[2] x[n];
}
parameters {
vector[2] mu;
vector<lower=0>[2] lambda;
real<lower=-1,upper=1> r;
}
transformed parameters {
vector<lower=0>[2] sigma;
cov_matrix[2] T;
Stan内の多変量正規分布は分散共分散行列
の代わりに共分散行列利用：多変量正規分布
は変数をベクトル化する必要がある。](https://image.slidesharecdn.com/rbayse2015710-150727085101-lva1-app6891/85/R-1-BCM-32-320.jpg)
![// Reparameterization
sigma[1] <- inv_sqrt(lambda[1]);
sigma[2] <- inv_sqrt(lambda[2]);
T[1,1] <- square(sigma[1]);
T[1,2] <- r * sigma[1] * sigma[2];
T[2,1] <- r * sigma[1] * sigma[2];
T[2,2] <- square(sigma[2]);
}
model {
// Priors
mu ~ normal(0, inv_sqrt(.001));
lambda ~ gamma(.001, .001);
// Data
x ~ multi_normal(mu, T);
}"
λの逆平方根＝σ
逆平方Gamma分布をモデルで
は仮定してた：逆平方根の処理
をしてからGamma関数仮定
※豆知識：Gamma分布
＝発生率1/βの事象が複数(α)
回起きるまでの待ち時間分布](https://image.slidesharecdn.com/rbayse2015710-150727085101-lva1-app6891/85/R-1-BCM-33-320.jpg)
搁でベイズをやってみよう！(コワい本1章)蔼叠颁惭勉强会
- 6. 実際の例①：BCM(コワイ本)の1章 ? 同じくらいの難しさの問題が10問あるテスト ? 知りたいことはあなたの能力 (θ＝正解率 ? 直接あなたの能力であるθは見れない 観察可能：テスト得点のみ
- 8. 例③ ? 実際のあなたの回答：10問中9問正解 ? 観察データ(D：9/10)によって → θ知識(一様分布なう)が更新 ? 事後分布p(θ│D)＝θに関する 各確率変数が真値の条件つき確率
- 9. ベイズの公式 ? べイズの公式 データ(＝尤度p(D│θ))+事前分布(＝p(θ)：一様分布) → 事後分布(＝ p(θ│D)) ? 言い換えれば事後＝尤度×事前/データが得られる確率 p(θ|D)=p(D|θ)p(θ)/p(D) 事後分布 事前分布
- 13. ベイズのいいところ：予測 ? 事後分布の使い方の一つ：予測 ? 前例と同じレベルの難しさの新しい問題5問 前回の事後分布(p(θ│n=10、k=9)を使って問題(nrep=5)をどう予測？ ? 数学的には事後を積分することで追加5問題の正解数(=krep)予測 ∫p(krep|θ,nrep=5)p(θ|n=10,k=9)dθ 事後分布で予測も可能！！
- 14. ベイズのいいところ２：逐次更新 ? 別個の情報をつなげたいときベイズは有用 例：同じレベルの難しさの問題5問 → 正解率：5問中3問 ? 過去に得た事後分布を事前分布に(θのアップデート方法) ? 観測するたびにデータが更新されていく(＝ベイズ更新)
- 15. ? 事後分布＝知りたいパラメータの不確実性記述 ＆ 確率的予測や逐次更新に有用 ? 尤度(データ)が二項分布に従う場合：一様事前＝ベータ 分布(α＝１、β＝１)とデータを組み合わせる→事後もま たベータ分布(α+k、β+ｎ-k) ? それらのような単一の結合例だと、事前と事後は同じ分 布になるけど、多くの場合そう都合よくならないよね。 → 複雑なデータには対応できない？ つまり上記の例は 単純すぎるって話
- 16. Markov Chain Monte Carlo ? （上の例の）事後分布は比較的単純データに対してしか使えない ＝ 現実の複雑なモデルはベイズ的解釈を使えなかった ? 上記の現状を打破したのがマルコフ連鎖モンテカルロ法と呼ばれる発 展的コンピュータ駆動のサンプリング手法(通称MCMC) → Gibbs サンプリング、Metropolis-Hastingsアルゴリズムのような MCMCを使い、事後分布を抽出可能に！ ? いまや“Bayesian models are limited only by the user’s imagination.”
- 17. MCMCって結局よくわかんない ? 本資料ではMCMCアルゴリズムについて詳しくは語りません(コワ イ本と同様に) ? とりあえずは一方向でない乱数シュミレーションを繰り返して妥 当な結果を出す道具と理解すれば大丈夫かと思います。 ? 詳しく知りたい人は https://www.youtube.com/watch?v=4gNpgSPal_8を参照した り、竹林先生の資料http://www.slideshare.net/yoshitaket/32- 35647139見たり、自分で検索してみてください。
- 18. ベイズのメリット?デメリットまとめ by清水先生 メリット ? データや推定値についての分布の制約がない(モデルが要請する制約 はある(例：回帰分析＝残差が正規分布 ? 事前分布を活用可能 デメリット ? 推定に時間がかかる ? 実際にどう使えばいいかわからない(→この会で分かる！はず)
- 19. Stan vs BUGS ?ここがすごいよStanくん！ ? No-U-Turn Sampler(NUTS)を利用：階層モデル、混合モデル、潜在 変数モデルなど変数間の相関が高くなる複雑なモデル向き ※BUGSはGibs Sampler, Metropolis-Hastingsアルゴリズム利用 ? さらにStanではBUGSでは許されない非正則な事前分布利用可能 ＋Shinystanなど関連する開発がいまもなお活発 心理データへの実用性が 高いということ
- 22. モデル ※[“”]はちゃんとコピペできないので書き直してね stancode<- “ data { int<lower=1> n; int<lower=0> k; } parameters { real<lower=0,upper=1> theta; } model { theta ~ beta(1, 1); // Observed Counts k ~ binomial(n, theta); }”
- 23. 続き library(rstan) datas <- list(n=10, k=9) # parameters to be monitored: parameters <- c("theta") # The following command calls Stan with specific options. samples <- stan(model_code=stancode, data=datas, pars=parameters, iter=20000, chains=2, thin=1, # warmup = 100, # Default = iter/2 # seed = 123 # Setting seed; Default is random seed )
- 24. Stan modelのすっごい適当な説明 stancode<- “ data { #ここのdataは実際に存在するデータを定義 int<lower=1> n; #問題の数nをつっこんでくれって意味：lowerは下限 int<lower=0> k; #正解数kを〃 } parameters { #知りたいパラメータ(データ上観察しえない)を定義 real<lower=0,upper=1> theta; #real=実数いれてやで＋<>内=下限～上限 → 能力θ } model { #分布の設定 theta ~ beta(1, 1); #事前分布：一様分布である(1，1)のβ分布 // Observed Counts k ~ binomial(n, theta); }” #正解数はパラメータ（n, θ）の二項分布に従うよ
- 25. Stan codeのすっごい適当な説明 samples <- stan(model_code=stancode, (さっき書いたモデル) data=datas, pars=parameters, iter=20000, #サンプリング数のこと。適切な数はデータによる chains=2, #ランダムサンプリングを行う回数 thin=1, #はしょる数：全部取り出すと多いのでiterを1つ飛びで抽出 # warmup = 100, #burn out期間：最初の方は収束してないのよね # Default = iter/2 #default # seed = 123 #乱数の指定; Default is random seed(＝これでいい) )
- 26. # The commands below are useful for a quick overview: print(samples) print(summary(samples)) 公式HPにはさらに結果 を見るためいろんなコ ードが実際にのってる けどそんな煩雑なこと をしなくてもShinystan というパッケージを使 えば結果がきれいに見れます (例)公式HPの全然関係ないデータ
- 29. 実際のデータにはどう使うの？ ? 上記のデータは尤度(データ)が二項分布に従う単純な例だった んだよね？ もっと慣れ親しんだもので やってみたい！ みんなご存じ相関係数を ベイズ推定してみよう！！
- 30. 相関係数を出そう！ ※BCM5章より ? ピアソンの積率相関係数 (r) ＝2変量の関係 －１～＋１ 計算式：r = xとyの共分散/SDx × SDy ? 通常は点推定値 (1値) として報告される r を ベイズ的に推定：相関の程度
- 31. Pearsonの積率相関係数のモデル ※詳細はCodeにて i data xi:観測データ 平均μ, 標準偏差σの多変量正規分布 r:相関係数 -1から1の範囲の無情報事前分布仮定 xi μ σ r μ1,μ2 ~ Gaussian(0,0.001) σ1,σ2 ~ InvSqrtGamma(0.001,0.001) r ~ Uniform(-1,1) xi ~ MvGaussian((μ1,μ2), σ2 1 rσ1σ2 ) rσ1σ2 σ2 2
- 32. Code：長いから覚悟しろよ！ ※ “” はコピペすると以下略 ? model <- " // Pearson Correlation data { int<lower=0> n; vector[2] x[n]; } parameters { vector[2] mu; vector<lower=0>[2] lambda; real<lower=-1,upper=1> r; } transformed parameters { vector<lower=0>[2] sigma; cov_matrix[2] T; Stan内の多変量正規分布は分散共分散行列 の代わりに共分散行列利用：多変量正規分布 は変数をベクトル化する必要がある。
- 33. // Reparameterization sigma[1] <- inv_sqrt(lambda[1]); sigma[2] <- inv_sqrt(lambda[2]); T[1,1] <- square(sigma[1]); T[1,2] <- r * sigma[1] * sigma[2]; T[2,1] <- r * sigma[1] * sigma[2]; T[2,2] <- square(sigma[2]); } model { // Priors mu ~ normal(0, inv_sqrt(.001)); lambda ~ gamma(.001, .001); // Data x ~ multi_normal(mu, T); }" λの逆平方根＝σ 逆平方Gamma分布をモデルで は仮定してた：逆平方根の処理 をしてからGamma関数仮定 ※豆知識：Gamma分布 ＝発生率1/βの事象が複数(α) 回起きるまでの待ち時間分布
- 34. Data x <- matrix(c( .8, 102, 1.0, 98, .5, 100, .9, 105, .7, 103, .4, 110, 1.2, 99, 1.4, 87, .6, 113, 1.1, 89, 1.3, 93), nrow=11, ncol=2, byrow=T) 行?列の指定
- 35. n <- nrow(x) # number of people/units measured data <- list(x=x, n=n) # to be passed on to Stan myinits <- list( list(r=0, mu=c(0, 0), lambda=c(1, 1))) # parameters to be monitored: parameters <- c("r", "mu", "sigma") # The following command calls Stan with specific options. # For a detailed description type "?rstan". samples <- stan(model_code=model, data=data, init=myinits, # If not specified, gives random inits pars=parameters, iter=10000, chains=1, thin=1, # warmup = 100, # Stands for burn-in; Default = iter/2 # seed = 123 # Setting seed; Default is random seed ) Init＝ランダムサンプリングす る時の最初の値の決定
- 36. おめでとう！やっとゴールやね！ ? launch_shinystan(samples) ? こんなんできるよ → だいたい-0.8くらいですね。 さらに事前確率やエラーバーを考慮 したより自由な分析も可能です。 ※詳しくはBCMの5章後半に
- 39. Bayes統計と従来の統計の違い ? Bayesは正規分布以外の分布を複数仮定し、事前情報を自由に 更新することができるという利点がある。 ? 一方で本当に自由な分析を可能にするためにはデータ?分布? 関数などの理解が必要不可欠である。 ? 従来の統計は、SPSS, Amosといった分析用のソフトでお手軽 にすることが可能であるということが利点の一つといえるだろ う（一方で統計解析がブラックボックスになりうる危険性も）。
- 41. おしまいです。 今日からベイズ！