Question

I am trying to find HPD credible set using R --

Sample values for X: 1, 2, 0, 1

Poisson Distribution, posterior is a gamma distribution with alpha = 0.5 and beta = 0.4

What is 95% HPD credible set using R?

Answer 1

Tried this concept for the first time. Let me know if it helped.

R-code

===============================================================================================

x <- c(1, 2, 0, 1)
lambda_true <- var(x)
n <- length(x)
alpha <- 0.5
beta <- 0.5
alpha_conf <- 0.05
# A posterior distribution for the parameter λ is Gamma ( ny'+α, n+β)
alpha_1 <- sum(x) + alpha
beta_1 <- n + beta

#compute 0.025- and 0.975-quantiles
q_lower <- qgamma(alpha_conf / 2, alpha_1, beta_1)
q_upper <- qgamma(1 - alpha_conf / 2, alpha_1, beta_1)
c(q_lower, q_upper)
# [1] 0.3000433 2.1136409 upper and lower bound)
# set up grid for plotting

# set up grid for plotting

lambda <- seq(0,7, by = 0.001)
lambda_true

plot(lambda, dgamma(lambda, alpha_1, beta_1), type = 'l', lwd = 2, col = 'violet',
     ylim = c(0, 1.5), xlab = expression(lambda),
     ylab = expression(paste('p(', lambda, '|y)')))

y_val <- dgamma(lambda, alpha_1, beta_1)
x_coord <- c(q_lower, lambda[lambda >= q_lower & lambda <= q_upper], q_upper)
y_coord <- c(0, y_val[lambda >= q_lower & lambda <= q_upper], 0)
polygon(x_coord, y_coord, col = 'pink', lwd = 2, border = 'violet')
abline(v = lambda_true, lty = 2)

lines(lambda, dgamma(lambda, alpha, beta),
      type = 'l', lwd = 2, col = 'orange')
legend('topright', inset = .02, legend = c('prior', 'posterior'),
       col = c('orange', 'violet'), lwd = 2)

===============================================================================================

0.3000433 2.1136409 upper and lower bound respectively,