# # Chapter 4 -- Bayesian Computation With R # Sleeping Habits Problem # # Remove all objects just to be safe # rm(list=ls(all=TRUE)) # library(LearnBayes) # data <- c(9.0, 8.5, 7.0, 8.5, 6.0, 12.5, 6.0, 9.0, 8.5, 7.5, 8.0, 6.0, 9.0, 8.0, 6.0, 7.0, 10.0, 9.0, 7.5, 5.0, 6.5) # # S is the sum of squares # S <- sum((data - mean(data))^2) n <- length(data) # # sigma2 is a 1000 random draw with entries equal to S divided by draws from the # chi-square distribution with n-1 degrees of freedom # sigma2 <- S/rchisq(1000, n - 1) # # mu is a 1000 random draw from a N(mean, sd) distribution # mu <- rnorm(1000, mean = mean(data), sd = sqrt(sigma2)/sqrt(n)) muquantile <- quantile(mu, c(0.05, 0.95)) sigmaquantile <- quantile(sqrt(sigma2), c(0.05, 0.95)) # prob75 <- mu + 0.647*sqrt(sigma2) meanprob75 <- mean(prob75) sdprob75 <- sd(prob75) # # Alternatively, use the normpostsim function in LearnBayes package # result <- normpostsim(data, m=1000) muquantilea <- quantile(result$mu, c(0.05, 0.95)) sigmaquantilea <- quantile(sqrt(result$sigma2), c(0.05, 0.95)) # prob75a <- mu + 0.647*sqrt(result$sigma2) meanprob75a <- mean(prob75) sdprob75a <- sd(prob75)