Нахождение задних интервалов для нескольких образцов

Question

Мне нужно сгенерировать 100 выборок размера 10 в R из гамма-распределения (10, 1), а затем вычислить 95% апостериорный интервал для каждого из них и построить их график. Последняя часть, с которой я борюсь. До сих пор я использовал:

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)

Но это не похоже на работу. Заранее спасибо, я здесь действительно не знаю.

Answer 1

Немного длинное решение, но у меня не было времени подумать над тем, чтобы сделать его более эффективным, поскольку я на работе :) Для меня это имеет смысл, если я не понял вашего вопроса ... Также для сюжета, вероятно, естьэто более чистый способ сделать это с помощью ggplot2, но я сделал это в цикле ... Последнее замечание: гамма-распределение имеет несколько различных параметризаций, которые вы можете использовать.Я предположил, что вы указываете Shape = 10, Rate = 1. Если есть вопросы, дайте мне знать!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")
}