library(LearnBayes)

###3.2関係
#data(footballscores) #footballscoresの読み出し
#attach(footballscores)
#d <- favorite - underdog - spread #得失点差の計算
#n <- length(d) #データ数の計算
#v <- sum(d^2) #総平方和の計算
#
#P <- rchisq(1000, n)/v #自由度nのΧ2分布からデータを1000個ランダムサンプリング。
#s <- sqrt(1/P)
#hist(s, main="") #これでp.41の図が出力

###3.3関係
#病院1の場合
alpha <- 16; beta <- 15174
yobs <- 1; ex <- 66
y <- 1: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)

#病院2の場合
ex <- 1767; yobs <- 4
y <- 0:10
py <- dpois(y, lam*ex)*dgamma(lam, shape = alpha, rate = beta)/dgamma(lam, shape = alpha + y, rate = beta + ex)
cbind(y, round(py, 3))
lambdaB <- rgamma(1000, shape = alpha + yobs, rate = beta + ex)

##病院1と2の結果の違いを図示
#図示する際のxの値の決定
lambda <- seq(0, max(c(lambdaA, lambdaB)), length=500)
#グラフィックウィンドウに描く図の配置を指定
par(mfrow=c(2, 1))
hist(lambdaA, freq = FALSE, main="", ylim=c(0, 1500))
lines(lambda, dgamma(lambda, shape = alpha, rate = beta))
hist(lambdaB, freq = FALSE, main="", ylim=c(0, 1500))
lines(lambda, dgamma(lambda, shape = alpha, rate = beta))


###3.4関係
##事前分布が正規分布の場合の計算
mu <- 100 #平均値（事前分布）
tau <- 12.16 #標準偏差（事前分布）
sigma <- 15 #各テストの点の標準偏差
n <- 4 #テスト回数
se <- sigma/sqrt(4) #テストの平均値の標準偏差（上記はすでに与えられているもの）
ybar <- c(110, 125, 140) #テストの点数
tau1 <- 1/sqrt(1/se^2 + 1/tau^2) #thetaの事後分布の標準偏差
mu1 <- (ybar/se^2 + mu/tau^2) * tau1^2 #thetaの事後分布の平均値
summ1 <- cbind(ybar, mu1, tau1) #各テストの点数、推定された平均値、標準偏差

##事前分布がt分布の場合に、正規分布のばらつきと対応させるためのスケールパラメータの計算
tscale <- 20/qt(0.95, 2) #自由度2のt分布における95%分位点

##事前分布としての正規分布とt分布を描画
theta <- seq(60, 140, length=200)
plot(theta, 1/tscale*dt((theta - mu)/tscale,2), type="l", ylab="Prior Density")
lines(theta, 1/10*dnorm((theta - mu)/tau), lwd=3)

##事前分布がt分布の場合の計算
summ2 <- c()
for (i in 1:3) {
     theta <- seq(60, 180, length=500) #平均値の範囲
     like <- dnorm((theta - ybar[i])/7.5) #尤度の計算
     prior <- dt((theta - mu)/tscale, 2) #事前分布の計算
     post <- prior * like #上記2つをかける
     post <- post/sum(post) #postの合計でpostを割ることで無理やり確率密度関数化
     m <- sum(theta * post) #平均値の算出
     s <- sqrt(sum(theta^2 * post) - m^2) #標準偏差の算出
     summ2 <- rbind(summ2, c(ybar[i], m, s)) #テストの点数、平均値、標準偏差を表にする
     }

##事前分布の違いによる事後分布の違いの作図
normpost <- dnorm(theta, mu1[3], tau1)
normpost <- normpost/sum(normpost)
plot(theta, normpost, type="l", lwd=3, ylab="Posterior Density")
lines(theta, post)
legend(locator(1), legend=c("t prior", "normal prior"), lwd=c(1,3))


###3.5関係
##
n <- 20
y <- 5
a <- 10
p <- 0.5
m1 <- dbinom(y, n, p) * dbeta(p, a, a)/dbeta(p, a + y, a + n - y)
lambda <- dbinom(y, n, p)/(dbinom(y, n, p) + m1)

##事前分布beta(a, a)のaを変化させた場合の作図
loga <- seq(-4, 5, length=100)
a <- exp(loga)
m2 <- dbinom(y, n, p) * dbeta(p, a, a)/dbeta(p, a + y, a + n - y)
lambda <- dbinom(y, n, p)/(dbinom(y, n, p) + m2)
plot(loga, lambda, type="l", xlab="log(a)", ylab="Prob(coin is fair)")

#観測された表の回数が5回未満の場合
n <- 20
y <- 5
a <- 10
p <- 0.5
m2 <- 0
for (k in 0:y) { #0:yによって表の回数が0回から5回までの確率を全て算出する
     m2 <- m2 + dbinom(k, n, p)*dbeta(p, a, a)/dbeta(p, a + k, a + n -k)
     }
lambda <- pbinom(y, n, p)/(pbinom(y, n, p) + m2)