# Coverage probability for the Bayesian credible interval for Normal distribution

Bayesian Inference for the Normal Distribution, I use the following r code to obtain the posterior distribution. Let's say the data, $$X\sim N(\mu, \sigma^{2})$$ and $$\mu \sim N(0,10)$$ and $$\sigma \sim \text{discrete uniform} (0,10)$$. I used rethinking package in r to build the above model. Also, I simulated the $$10000$$ samples from the posterior distribution. It can be done as follows.

install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')
library(rethinking)

set.seed(6)
x <- data.frame(x = rt(100,3))
xbar <- mean(x\$x)
fit <- rethinking::map(
alist(
x ~ dnorm(mu, sigma),
mu ~ dnorm(1, 10),
sigma ~ dunif(0, 50)
),
data=x)

precis(fit, corr=TRUE)

sim_post <- extract.samples(fit)
dim(sim_post)
post_par <- apply(sim_post, 2, mean)
post_par

quantile(sim_post$$mu , c(.05, .95)) quantile(sim_post$$sigma, c(0.05, 0.95))

level <- 0.95

CI_low <- qnorm((1-level)/2, sim_post[,1], sim_post[,2])
CI_up <- qnorm((1+level)/2, sim_post[,1], sim_post[,2])

cover_prob <- mean((CI_low < 0) & (0 < CI_up))
cover_prob


I graphed both posterior parameters below.

How can I obtain the coverage probability for mu? Why am I getting $$1$$ always? Are there any problem in r code? Are there any ways to obtain the coverage probability for this analysis?

• It's not clear what kind of object sim_post might be and therefore makes it difficult to determine how you are trying to compute the interval. Nevertheless, I believe the answers to your questions will immediately become apparent if you were to visualize the posterior distribution: why not plot its histogram, for instance? – whuber Sep 26 '18 at 20:49
• @whuber, the dim(sim_post) is [1] 10000 2, it contains simulated values of the posterior parameters. – score324 Sep 26 '18 at 20:51