How does one calculate a 95% Credible Interval for Bernouli sampling? I have looked into the package emdbook and the function tcredint() but i am unsure if this the way to go.

If anyone has references, I'm grateful.

Create data with 2000 rows of (0,1) set.seed(100) data<-data.frame(replicate(1,sample(0:1,2000,replace=TRUE))) colnames(data) <- c("sample") require(psych) mean(data$sample) sd(data$sample)'

Confidence Interval require(Publish) ci.mean(data)]

Credible Interval install.packages("gmodels") require(gmodels) require(stats) p=pbeta(data$sample, shape1=2, shape2=2, ncp = 0, lower.tail = TRUE, log.p = FALSE) ci(p)


1 Answer 1


Technically, you've not supplied enough information to answer as we'd need to know your prior.

However, I will assume you're looking at a $p\sim$Beta$(\alpha, \alpha)$ prior. In that case, if you have $N$ trials with a total of $Y$ successes the posterior distribution is $p | \mathcal D \sim$Beta$(\alpha + N - Y, \alpha + Y)$.

You can use the distribution function for Beta to calculate the necessary tails using those parameters.

  • 1
    $\begingroup$ My prior is πθ~ß(a,b). From what I am working with, it looks like I am using a Beta-Binomial posterior distribution (π(θ/y) ~ ß (y+a,n–y+b)). So, my best bet is the pbeta function in R? $\endgroup$
    – Bryan Mac
    May 21, 2018 at 15:38
  • 1
    $\begingroup$ Yup. One detail 'Beta-Binomial' also refers to the predictive distribution for a new set of trials. You're just intereseted in the Beta posterior on the parameter. $\endgroup$ May 21, 2018 at 15:58
  • $\begingroup$ When i used the function pbeta, the credible intervals' lower and upper bound is identical to the confidence interval. Is my R script above accurate? I'm assuming that the reason is because of my Shape1 and Shape2 parameters. $\endgroup$
    – Bryan Mac
    May 22, 2018 at 2:57
  • $\begingroup$ Are you adding in the prior alpha and beta to the conjugate beta posterior shape parameters? $\endgroup$ May 22, 2018 at 22:37
  • $\begingroup$ I don't think I did unless the function, p=pbeta(data$sample, shape1=2, shape2=2, ncp = 0, lower.tail = TRUE, log.p = FALSE), did it for me. Is there an R script that I would need to do before calculating the pbeta? $\endgroup$
    – Bryan Mac
    May 23, 2018 at 2:50

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.