Finding posterior intervals for multiple samples

Question

I have to generate 100 samples of size 10 in R from a Gamma (10, 1) distribution and then compute the 95% posterior interval for each of them and plot them. The last part I struggle with. So far I used:

data = rgamma(n,10,1)
simdata = rgamma(n=100, 10, 1)
matrixdata = matrix(simdata, nrow=10, ncol=10)
quantile(matrixdata, HPD=TRUE, MM=TRUE, prob=0.95, plot=FALSE, PDF=FALSE)

But that doesn't seem to work. Thank you in advance, I'm really clueless here.

Answer 1

Bit of a lengthy solution but did not have the time to consider making it more efficient as I am at work :) To me, this makes sense, unless I misunderstood your question... Also for the plot, there is probably a way more clean way to make it using ggplot2, but I did it in a loop... Final note, the Gamma distribution has several different parametrizations you can use. I assumed you specify Shape = 10, Rate = 1. If any questions let me know! Peaaace

no_simulations <- 100
n <- 10
shape <- 10
rate <- 1
set.seed(10)
empirical_quantile_intervals <- matrix(ncol = 2, nrow = no_simulations)
names(empirical_quantile_intervals) <- c("Q025", "Q975")

simulation_matrix <- matrix(nrow = no_simulations, ncol = n)
for (i in 1:no_simulations) {
  simulation_matrix[i, ] <- rgamma(n = n, shape = shape, rate = rate)
  empirical_quantile_intervals[i, 1] <- quantile(simulation_matrix[i, ], probs = 0.025)
  empirical_quantile_intervals[i, 2] <-  quantile(simulation_matrix[i, ], probs = 0.975)
}

plot(empirical_quantile_intervals[1, ], c(1,1), xlim = c(0, max(empirical_quantile_intervals)), 
     ylim = c(0, no_simulations), type="b", 
     xlab = "Quantile intervals", 
     ylab = "Simulation")

for(i in 2:no_simulations) {
  lines(empirical_quantile_intervals[i, ], c(i,i), type="b")
}