# calculating credible intervals for Bernoulli trials with a beta-distribution

I'm interested in finding credible intervals for Bernoulli trials with a variable probability of success.

So, the process (as I understand it) is the probability of success, p, is modeled with a beta distribution based on an expected outcome and your experimental results:

$$p ∼ beta(\alpha, \beta)$$

This value of p is then feed into

$$x ∼ bernoulli(p)$$

and from this credible intervals can be generated. Actually doing this calculation and finding the credible intervals is where I'm having trouble.

I've seen a few different ways of doing this

I have some results from some preliminary experiments I've done; $$\alpha$$ = 51 + 1, $$\beta$$ = 9 + 1, n = 60 (with $$\alpha_0 = \beta_0 = 1$$).

Link 1 method: (0.82, 0.85). The bootstrap results are essentially the same. This interval seems very small.

Link 2 method: (0.74, 0.92) which is calculated by qbeta(c(0.025, 0.975), 51 + 1, 9 + 1) Seems to ignore using p to feed back into x ~ bernoulli(p), and just calculates credible intervals directly from the beta distribution.

Link 3 method: I conceptually understand but am having a hard time putting it into code.

Should the interval for Link 1 be as small as it is? Is qbeta skipping x ~ bernoulli(p), and does that matter? Any help with calculating this would be much appreciated.

This is the first time I've used R. I have some experience using Python.