ベイズ統計 第2回演習

  • ここでは、慶應大学SFCでの開講科目「ベイズ統計」での演習で用いている分析例やRコード・JAGSコードなどを、授業履修者のために示しています。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
  • Rパッケージやデータの更新などにより、本書執筆時とはデータやRコードを予告無しに変更する場合があります。 本書記載のR分析例以外は、R最新版等で動作確認できていないRコードも存在しますが、参考のために記載しておきます(旧バージョンで動作することがあります)。

第2回R演習の内容

  • ベイズ更新と経験ベイズのシミュレーションについて理解を深めます。
  • 授業第6回・第7回の講義資料を適宜参照してください。

自然共役事前分布 二項分布とベータ分布

  • 自然共役事前分布が二項分布のとき、尤度にベータ分布を選ぶことで、ベータ分布に従う事後分布が得られます。
  • 過去に似たようなチームと対戦して10戦6勝(n=10, r=0.6)だったとすると、その二項分布は以下のようになる。
    # Βinom(10,0.6)
    x <- seq(0,20,by=1)
    plot(x, dbinom(x, size=10, prob=0.6), type="l", lwd=2)
  • Βinom(10,0.6)はΒ(7,5) (α=7, Β=5) に比例する。⇒尤度となる
    x <- seq(0,1,by=0.01)
    par(family="HiraMaruProN-W4")
    plot(x, dbeta(x, 7, 5), type="l", ylim=c(0,4), lwd=2, xlab="密度", ylab="確率", col="grey")
  • 尤度が最大となる値を得ると、r=0.6のときに尤度=2.759となることがわかる。
    # 数値計算の解から見つけ出す方法
    lf <- cbind(x,dbeta(x,7,5))
    lf[lf[,2]==max(dbeta(x,7,5))]
    ## 尤度関数を最大化する方法
    # L <- function(x){x^6*(1-x)^4} # 二項分布を尤度関数に採用
    # optimize(L, c(0,1), maximum=TRUE)
  • 事前分布は、理由不十分の原則より、B(2,2)
    curve(dbeta(x, 2, 2), type="l" , col="red", add=TRUE, lwd=2)
    # B(2,2)は0.5で最大となる
  • 事後分布はΒ(8,6) (α=8, Β=6) に比例する。
    curve(dbeta(x, 8, 6), type="l" , col="blue", lwd=2, add=TRUE)
  • 凡例の表示
    legend(.7,4,c("尤度関数","事前分布","事後分布"), col=c("grey","red","blue"),lwd=c(2,2,2))

自然共役事前分布◆Д櫂▲愁麒布とガンマ分布

  • 講義では、事前分布と尤度がポアソン分布に従う例として、サッカーJリーグの各チームの一試合辺り得点分布を挙げました。
  • プレシーズンマッチ5試合の結果のように、事前分布が平均1.38、分散0.85のポアソン分布に従うとき、ガンマ分布Γ(2.24, 162)となる(α=2.24、λ=1.62)。
    x <- seq(0,10,by=0.01)
    par(family="HiraMaruProN-W4")
    plot(x, dgamma(x, shape=2.24, scale=1.62), type="l", lwd=2, ylim=c(0,0.6), xlab="得点", ylab="確率")
  • 一方、2013年のJリーグ一チーム一試合あたりの得点は平均1.44、分散1.38のポアソン分布に従う事がわかっており、これをガンマ分布で表すとΓ(1.50, 1.04)となります(α=1.50、λ=1.04)。
    curve(dgamma(x, shape=1.50, scale=1.04), type="l", lwd=2, col="red", add=TRUE)
  • すると事後分布は、Γ(2.74, 2.66)となる。
    curve(dgamma(x, shape=2.74, scale=2.66), type="l", lwd=2, col="blue", add=TRUE)
    legend(7,.5,c("事前分布","尤度関数","事後分布"), col=c("black","red","blue"),lwd=c(2,2,2))

自然共役事前分布:正規分布のベイズ推定(平均未知・分散既知)

  • アメリカンフットボールNFLの2011-2013年シーズンのデータから、ホームチームの得点とアウェイチームの得点の差が正規分布に従うと考えられる。実際の得点分布などのデータは、第1回R演習からダウンロード可能です。
  • 2011-2013年の3年間合計の試合別得点差分布はN(5.42, 13.54)に従うと考えられる。
    x <- seq(-40,60,by=0.01)
    par(family="HiraMaruProN-W4")
    # 確率分布
    plot(x, dnorm(x, mean=5.42, sd=13.54), type="l", lwd=2, ylim=c(0,0.05), xlab="得点差", ylab="確率")
  • 2011年、2012年、2013年の試合別得点差分布はそれぞれ、N(5.66, 13.74)、N(4.91, 14.44)、N(5.70, 12.40)に従うと考えられる。
    curve(dnorm(x, mean=5.66, sd=13.74), type="l", lwd=2,col="purple", add=TRUE)
    curve(dnorm(x, mean=4.91, sd=14.44), type="l", lwd=2,col="lightgreen", add=TRUE)
    curve(dnorm(x, mean=5.70, sd=12.40), type="l", lwd=2,col="orange", add=TRUE)
    legend(25,.04,c("2011-2013年","2011年","2012年","2013年"), col=c("black","purple","lightgreen","orange"),lwd=c(2,2,2,2))
  • 事後正規分布(平均)の推定
    # 事前分布の平均と標準偏差
    mu0 <- 6.00 # 平均
    sig0 <- 13.62 # 標準偏差
    # 尤度関数の平均と標準偏差及びサンプル数
    xvar <- 5.42 # 平均
    sig <-  13.54 # 標準偏差
    n <- 3 # サンプル数
    # 事後分布の平均と標準偏差及び精度
    mu1 <- ((mu0/sig0^2)+(n*xvar/sig^2))/((1/sig0^2)+(n/sig^2))
    tau <- (1/sig0^2)+(n/sig^2)
    sig1 <- sqrt(1/tau)
    # 事後正規分布
    x <- seq(-40,60,by=0.01)
    par(family="HiraMaruProN-W4")
    plot(x, dnorm(x, mean=mu1, sd=sig1), type="l", lwd=2, ylim=c(0,0.06), xlab="得点差", ylab="確率")
    # 尤度関数
     curve(dnorm(x, mean=xvar, sd=sig), type="l", lwd=2,col="red", add=TRUE)
    # 事前分布
     curve(dnorm(x, mean=mu0, sd=sig0), type="l", lwd=2,col="grey", add=TRUE)
    # 凡例
    legend(30, 0.06, c("事前分布", "尤度関数", "事後分布"), col=c("grey","red","black"), lwd=c(2,2,2))
  • 正規分布の平均をベイズ更新により推定する
    # 事前分布の平均と標準偏差
    mu0 <- 6.00 # 平均
    sig0 <- 13.62 # 標準偏差
    # 尤度関数の平均と標準偏差及びサンプル数
    xvar <- 5.42 # 平均
    sig <-  13.54 # 標準偏差
    n <- 3 # サンプル数
    # ベイズ更新の繰り返し回数
    n.update <- 100
    # ベイズ更新の結果を格納
    post.mu <- c(xvar,mu1)
    for(i in 1:n.update){
    	mu1 <- ((mu0/sig0^2)+(n*xvar/sig^2))/((1/sig0^2)+(n/sig^2))
    	tau <- (1/sig0^2)+(n/sig^2)
    	sig1 <- sqrt(1/tau)
    	# ベイズ更新の結果を格納
    	if(i==1){post.mu<-mu1}
    	if(i>1){post.mu <- c(post.mu,mu1)}
    	# print(mu1)
    	mu0 <- mu1
    	sig0 <- sig1
    	}
    # 事後平均を表示
    plot(post.mu, type="b")
    # 最後の10回分を表示
    tail(post.mu, 10)

自然共役事前分布ぁЮ亀分布のベイズ推定(平均未知・分散未知)

  • 事後正規分布(平均と分散)の推定
    # 事前分布の平均と標準偏差
    mu0 <- 6.00 # 平均
    sig0 <- 13.62 # 標準偏差
    # 尤度関数の平均と標準偏差及びサンプル数
    xvar <- 5.42 # 平均
    tau <- rgamma(1, shape=1, rate=1)
    sig <- sqrt(1/tau)
    n <- 3 # サンプル数
    # 事後分布の平均と標準偏差及び精度
    mu1 <- ((mu0/sig0^2)+(n*xvar/sig^2))/((1/sig0^2)+(n/sig^2))
    tau <- (1/sig0^2)+(n/sig^2)
    sig1 <- sqrt(1/tau)
  • ベイズ更新による事後平均と事後分散(標準偏差)の推定
    # 事前分布の平均と標準偏差
    mu0 <- 6.00 # 平均
    sig0 <- 13.62 # 標準偏差
    # 尤度関数の平均と標準偏差及びサンプル数
    xvar <- 5.42 # 平均
    tau <- rgamma(1, shape=1, rate=1)
    sig <- sqrt(1/tau)
    n <- 3 # サンプル数
    # ベイズ更新の繰り返し回数
    n.update <- 100
    for(i in 1:n.update){
    	mu1 <- ((mu0/sig0^2)+(n*xvar/sig^2))/((1/sig0^2)+(n/sig^2))
    	tau <- (1/sig0^2)+(n/sig^2)
    	sig1 <- sqrt(1/tau)
    	if(i==1){post.mu<-mu1; post.sig<-sig1}
    	if(i>1){post.mu <- c(post.mu,mu1)
    			post.sig <- c(post.sig,sig1)}
    	# print(mu1)
    	mu0 <- mu1
    	sig0 <- sig1
    	}
    # 事後平均を表示
    par(family="HiraMaruProN-W4", mfrow=c(2,1))
    plot(post.mu, type="b", xlab="更新回数", ylab="事後平均")
    plot(post.sig, type="b", xlab="更新回数", ylab="事後標準偏差")
    # 最後の10回分を表示
    tail(post.mu, 10)
    tail(post.sig, 10)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-06-23 (月) 10:05:50 (1944d)