Bayesian user survey with a credible interval [R code]

bayesian.R

> rdirichlet <- function(a) {
+     y <- rgamma(length(a), a, 1)
+     return(y / sum(y))
+ }
> x <- c(65, 71, 532, 307, 369, 234, 584)
> SIMS <- 10000
> t <- matrix(nrow = SIMS, ncol = length(x), data = 0)
> for (i in 1:SIMS) t[i,] = rdirichlet(x + 1)
> sum(t[,3]) / SIMS # estimate of "college or grad"
[1] 0.245605
> quantile(t[,3], probs = c(0.025, 0.975)) # credible interval for "college or grad"
     2.5%     97.5% 
0.2279062 0.2640242 
> sum(t[,3] > t[,7]) / SIMS # more "college or grad" than "other"
[1] 0.0588
> sum(t[,1] > t[,2]) / SIMS # more "k-12 teacher or librarian" than "k-12 student"
[1] 0.3076

# http://stats.stackexchange.com/questions/39319/bayesian-user-survey-with-a-credible-interval#40403