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

  1. https://stats.stackexchange.com/a/82640/197669

  2. https://rpubs.com/sitaramgautam/145048

  3. https://stats.stackexchange.com/a/238582/197669

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.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.