ベイズ統計R演習2

注意事項

• 慶應義塾大学SFCで開講している「ベイズ統計」の授業履修者向けの演習用ページです。
• 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
• Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

二項分布とベータ分布

• 自然共役事前分布が二項分布のとき、尤度にベータ分布を選ぶことで、ベータ分布に従う事後分布が得られます。
• 事前分布として、過去の対戦結果を採用し、１０戦６勝（n=10, r=0.6）だったとすると、その二項分布は以下のようになる。
• 二項分布とベータ分布の確率密度関数は、それぞれdbinom()とdbeta()

# Β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")
# 事前分布は、理由不十分の原則より、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))

• 尤度が最大となる値を得ると、r=0.6のときに尤度=2.759となることがわかる

# 数値計算の解から見つけ出す方法
lf <- cbind(x,dbeta(x,7,5))
lf[lf[,2]==max(dbeta(x,7,5))]

## [1] 0.600000 2.759049

# 尤度関数を最大化する方法
L <- function(x){x^6*(1-x)^4} # 二項分布を尤度関数に採用
optimize(L, c(0,1), maximum=TRUE)

## $maximum
## [1] 0.600002
## 
## $objective
## [1] 0.001194394

ポアソン分布とガンマ分布

• 講義では、事前分布と尤度がポアソン分布に従う例として、サッカーJリーグの各チームの一試合辺り得点分布を挙げました
• プレシーズンマッチ５試合の結果のように、事前分布が平均1.38、分散0.85のポアソン分布に従うとき、ガンマ分布Ga(2.24, 162)となる（α=2.24、λ=1.62）。
• 一方、2013年のJリーグ一チーム一試合あたりの得点は平均1.44、分散1.38のポアソン分布に従う事がわかっており、これをガンマ分布で表すとGa(1.50, 1.04)となります（α=1.50、λ=1.04）
• すると事後分布は、Ga(2.74, 2.66)となる

# 事前分布を描画
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="確率")
# 尤度を描画
curve(dgamma(x, shape=1.50, scale=1.04), type="l", lwd=2, col="red", add=TRUE)
# 事後分布を描画
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年シーズンのデータから、ホームチームの得点とアウェイチームの得点の差が正規分布に従うと考えられる。
• 2011-2013年の3年間合計の試合別得点差分布はN(5.42, 13.54)に従うと考えられる。
• 2011年、2012年、2013年の試合別得点差分布はそれぞれ、N(5.66, 13.74)、N(4.91, 14.44)、N(5.70, 12.40)に従うと考えられる。
• それぞれの分布は以下のように表せる。

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="確率")
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)

##  [1] 5.422092 5.422069 5.422047 5.422026 5.422004 5.421983 5.421963
##  [8] 5.421943 5.421924 5.421904

正規分布のベイズ推定（平均未知・分散未知）

• 事後正規分布（平均と分散）の推定

# 事前分布の平均と標準偏差
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))
mu1

## [1] 5.421619

# 事後分布の精度
tau <- (1/sig0^2)+(n/sig^2)
tau

## [1] 1.93148

# 事後分布の標準偏差（事後標準偏差）
sig1 <- sqrt(1/tau)
sig1

## [1] 0.7195398

• ベイズ更新による事後平均と事後分散（標準偏差）の推定

# 事前分布の平均と標準偏差
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)

##  [1] 5.420007 5.420007 5.420006 5.420006 5.420006 5.420006 5.420006
##  [8] 5.420006 5.420006 5.420006

# 事後標準偏差の最後の10回分を表示
tail(post.sig, 10)

##  [1] 0.04600758 0.04575686 0.04551019 0.04526747 0.04502859 0.04479346
##  [7] 0.04456197 0.04433403 0.04410956 0.04388846