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]