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.