Bayes Book
Chapter 5
R code
Adam Branscum





Section 5.1.1.2

# Code 1
theta0 <- 0.2
btilde <- 1:100
a <- (1+theta0*(btilde-2)) / (1-theta0)
qbeta(0.95,a,btilde)



# Code 2
btilde <- 1:100
theta0 <- 0.2
b1tilde <- 10+(btilde/100)
a <- (1+theta0*(b1tilde-2)) / (1-theta0)
qbeta(0.95,a,b1tilde)





# Section 5.2:  Brass alloy zinc data
n=12
y=c(4.20, 4.36, 4.11, 3.96, 5.63, 4.50, 5.64, 4.38,4.45, 3.67, 5.26, 4.66)
mean(y)
sd(y)
sd(y)/sqrt(n)
min(y)
max(y)
quantile(y, c(0.25,0.50,0.75))





# Section 5.2.3:  Independence priors
sigma0 <- 15.625
f <- 1:100
e <- 1+sigma0*f
qgamma(0.95,e,f)   ## Look for 20.31

f <- seq(2,3,0.001)
e <- 1+sigma0*f
qgamma(0.95,e,f)   ## Look for 20.31





tau0 <- 0.004096
d <- 1000:2000
c <- 1+tau0*d
qgamma(0.95,c,d)   ## Look for 0.008359




x <- seq(0,40,0.01)
y <- 1/sqrt(rgamma(10000,6.439,1328))
hist(y,prob=T)
lines(x,dgamma(x,41.375,2.588),lty=2)





# Exercise 5.26
#1 R Code for finding prior on sigma
alpha <- 0.90
beta <- 0.95
a <- 65 # Best guess for mu
tildegamma <- 85 # Best guess for gamma_alpha
tildeu <- 91 # Best guess percentile of gamma_alpha
zalpha <- 1.28  # qnorm(0.90,0,1)
f <- 3 # Initial value for f
# Could use a sequence of values, say f <- seq(1,50,1)
sigma0 <- (tildegamma - a)/zalpha
e <- 1 + sigma0*f
# We must find the Gamma(e,f) distribution that
# has beta-percentile = tildesigmabeta
tildesigmabeta <- (tildeu - a)/zalpha
trialq <- qgamma(beta,e,f) # Return beta-percentile for the selected gamma distribution.
trialq # If trialq = tildesigmabeta
tildesigmabeta # stop and pick corresponding f



#2 R Code for finding prior on tau
alpha <- 0.9
beta <- 0.95
zalpha <- 1.28
a <- 65
tildegamma <- 85
tildel <- 79
d <- 1100
tau0 = (zalpha/(tildegamma - a))^2
c = 1 + tau0*d
tildetaubeta = (zalpha/(tildel - a))^2
trialq = qgamma(beta,c,d)
trialq
tildetaubeta









#Exercise 5.32
mode <- 10
b <- 1:100
a <- 1+(mode * b)
qgamma(0.95,a,b)

mode <- 10
b <- 1:100/100
a <- 1+(mode*b)
round(qgamma(0.95,a,b),1)

b[22]
a[22]