Rでベイズをやってみる的なもの色々

来週のゼミが僕が担当なんだけど*1、数式展開して、「はい、こうなりますね」だけじゃ面白くないので、Rでやってみる。Bayesian Computation with R (Use R)をなぞっているだけですけどね。

正規分布で、平均既知、分散未知の場合

P30くらいからの内容。PRML読書会でもやったように、普通はガンマ分布を共役事前分布としておくと都合がよい。

ただ、Bayesian Computation with R (Use R)のほうでは、 $p(\sigma^2)=\frac{1}{\sigma^2}$ を事前分布に仮定している。いわゆる変則事前分布となっている*2。しかし、ƒKƒ“ƒ}[•ª•z‚Ì«Ž¿の付近とかを見てみると、ガンマ分布はパラメータをいじれば、この分布と見なせないこともないので、まあ問題ないとしておきましょう。精度パラメータ(分散の逆数)をPとすると、PはU/vに分布するらしい。Uは自由度nのX^2分布になるので(あとで確かめる)、じゃあ、乱数生成してヒストグラムを書けるよねという話になっていく。で、その事後分布を書いてみると以下のようになりますという感じ。

install.packages("LearnBayes")
library(LearnBayes)
library(help=LearnBayes)

data(footballscores)
footballscores

attach(footballscores)

d <- favorite - underdog - spread

n <- length(d)
v <- sum(d^2)

P <- rchisq(1000,n) / v
s <- sqrt(1/P)

hist(s)
quantile(s,probs=c(0.025,0.5,0.975))

指数分布

lambdaを推定しちゃうよ!!

alpha <- 16
beta <- 15174
yobs <- 0
ex <- 66

y <- 0:10
lam <- alpha / beta
py <- dpois(y,lam*ex) * dgamma(lam,shape=alpha,rate=beta) / dgamma(lam,shape=alpha+y,rate=beta+ex)
cbind(y,round(py,3))

lambdaA <- rgamma(1000,shape=alpha+yobs,rate=beta+ex)

hist(lambdaA,freq=FALSE)

正規分布で、平均、分散未知の場合

両方推定しちゃうよー。

data(marathontimes)
attach(marathontimes)
normchi2post

d <- mycontour(normchi2post,c(200,330,500,9000),time)
title(xlab="mean",ylab="variance")

S <- sum((time - mean(time))^2)

n <- length(time)
sigma2 <- S /rchisq(1000,n-1)
mu <- rnorm(1000,mean=mean(time),sd=sqrt(sigma2/sqrt(n)))
points(mu,sigma2)

quantile(mu,c(0.025,0.975))
quantile(sqrt(sigma2),c(0.025,0.975))