Bayesian analysis
R イ
難波 修史
イ ？ ！？
わ ！
資料 内容 本 第1章 内容
多分 含 了承
ふ わ イ 考え方 説明
：頻度論 従来 考え方 対比
従来 確率 頻度論
→無限回 試行 前提 確率
＝真 値 1 値 神 知
イ 確率
→ 時点 有 ー 確率 仮定
＝実際 ー ＋予備知識 確率
→確率的 真 値 1 値
イ 基本原理
．直観的信頼度 確率変数 定量化
．観察 ー 使 事前 情報 事後 情報
That’s it( )!
実際 例 ：BCM( ワイ本) 1章
同 難 問題 10問あ
知 あ 能力
(θ＝正解率
直接あ 能力 あ θ 見
観察可能： 得点
例
あ 能力θ 関 事前 情報(事前分布) 特定
※θ＝0～1 問題 関 情報(難 熟知度) 一切 状態
事前分布 p(θ)
※事前 情報 一切 (予測不可)
右 図 一様分布(全確率 同 )
例
実際 あ 回答：10問中9問正解
観察 ー (D：9/10) → θ知識(一様分布 ) 更新
事後分布p(θ│D)＝θ 関
各確率変数 真値 条件 確率
イ 公式
イ 公式
ー (＝尤度p(D│θ))+事前分布(＝p(θ)：一様分布)
→ 事後分布(＝ p(θ│D))
言 換え 事後＝尤度×事前/ ー 得 確率
p( θ|D)=p(D|θ)p(θ)/p(D)
事後分布 事前分布
事後 尤度×事前 比例関係 竹林先生 資料
Q ：
A ：事前分布 予備知識 実際
ー 事後分布
予測
例
F
igure1.1 参照
事前分布(Prior)
9/10 正解率 ー 情報
利用 事後分布(Posterior)
求 ：右図参照
新 ー 与え 分θ 不確実性 下 事前分
布 事後分布 狭 分布 ＝分布 尖
イ ：予測
事後分布 使 方 一 ：予測
前例 同 難 新 問題5問
前回 事後分布(p(θ│n=10 k=9) 使 問題(nrep=5) 予測？
数学的 事後 積分 追加5問題 正解数(=krep)予測
∫p(krep|θ,nrep=5)p(θ|n=10,k=9)dθ
事後分布 予測 可能！！
イ ：逐次更新
別個 情報 イ 有用
例：同 難 問題5問 → 正解率：5問中3問
過去 得 事後分布 事前分布 (θ アッ ー 方法)
観測 び ー 更新 ( ＝ イ 更新)
事後分布＝知 ー 不確実性記述
＆ 確率的予測 逐次更新 有用
尤度( ー ) 二項分布 従 場合：一様事前＝ ー
分布( ＝ ＝ ) ー 組 合わ →事後
ー 分布( +k +ｎ-k)
単一 結合例 事前 事後 同 分
布 多 場合 都合
→ 複雑 ー 対応 ？
単純 上記 例話
Markov Chain Monte Carlo
上 例 事後分布 比較的単純 ー 対 使え
＝ 現実 複雑 イ 的解釈 使え
上記 現状 打破 連鎖 ン カ 法 呼 発
展的 ン ュー 駆動 ン ン 手法(通称MCMC)
→ Gibbs ン ン Metropolis-Hastingsア MCMC 使 事後分布 抽出可能 ！
“Bayesian models are limited only by the user’s imagination.”
MCMC 結局 わ
本資料 MCMCア 詳 語 ( ワ
イ本 同様 )
あえ 一方向 乱数 ュ ー ョン 繰 返 妥
当 結果 出 道具 理解 大丈夫 思
詳 知 人
https://www.youtube.com/watch?v=4gNpgSPal_8 参照
竹林先生 資料http://www.slideshare.net/yoshitaket/32-
35647139見 自分 検索
イ ッ ・ ッ
by 清水先生
ッ
ー 推定値 分布 制約 ( 要請 制約
あ (例：回帰分析＝残差 正規分布
事前分布 活用可能
ッ
推定 時間
実際 使え わ (→ 会 分 ！ )
Stan vs BUGS
Stan ！
No-U-Turn Sampler(NUTS) 利用：階層 混合 潜在
変数 変数間 相関 高 複雑 向
※BUGS Gibs Sampler, Metropolis-Hastingsア 利用
Stan BUGS 許 非正則 事前分布利用可能
＋Shinystan 関連 開発 活発
心理 ー 実用性
高
Rstan 導入
ッ ー CRAN ウン ー
出来 (多分)
3行 ー 走
( C＋環境 同時 用意 必要有)
導入 詳細 小杉先生 イ
http://www.slideshare.net/KojiKosugi/r-stan あ 本家 イ (Rstan 検索検索 )
参照 難波 相談 意味
R 例
実際 ！！
※[ “”] 書 直
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); }”
続
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 )
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, θ 二項分布 従
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(＝ ) )
# The commands below are useful for a quick overview: print(samples)
print(summary(samples))
公式HP 結果
見
ー 実際
煩雑
Shinystan
ッ ー 使
え 結果 見 (例)公式HP 全然関係 ー
見 ！Shinystan( 導入)
https://github.com/stan-dev/shinystan/wiki
イン ー
devtoolsイン ー (多分 CRAN 思
)
→ library(devtools)
source_url("https://github.com/stan-
dev/shinystan/raw/develop/install_shinystan.R") install_shinyStan()
詳 上記 本家 イ 参照
結果 見方
library(shinyStan)
Sys.setlocale(locale=“English”) launch_shinystan(samples)
ｋ ※“” 注意
自分 書 直 終わ
Sys.setlocale(locale=“Japanese”)
実際 ー 使 ？
上記 ー 尤度( ー ) 二項分布 従 単純 例
？
慣 親 ！
存 相関係数
イ 推定 ！！
相関係数 出 ！
※BCM5 章
アソン 積率相関係数 (r) ＝2変量 関係
－ ～＋
計算式：r = x y 共分散/SDx × SDy
通常 点推定値 (1値) 報告 r イ 的 推定：相関 程度
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), σ21 rσ1σ2 ) rσ1σ2 σ22
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内 多変量正規分布 分散共分散行列
代わ 共分散行列利用：多変量正規分布
変数 化 必要 あ
// 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/ 事象 複数( )
回起 待 時間分布
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)
行・列 指定
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＝ ン ン ン
時 最初 値 決定
！ ー ！
launch_shinystan(samples)
→
-0.8
事前確率 エ ー ー 考慮 自由 分析 可能
※詳 BCM 5章後半
普通 Cor.test 1 行
Bayes 統計 従来 統計 違
Bayes 正規分布以外 分布 複数仮定 事前情報 自由
更新 利点 あ
一方 本当 自由 分析 可能 ー ・分布・
関数 理解 必要不可欠 あ
従来 統計 SPSS, Amos 分析用 ソ 手軽
可能 あ 利点 一 え
一方 統計解析 ッ ボッ 危険性
・・・ イ 利用 User次第
ベイ ズ
今日 イ ！
・・・
長 イ 的 相関 求 方
function 定義 使 思
