Suppose I flip a coin 584,369 times and observe 5,257 heads and 579,112 tails. I wish to know the lower edge of the 90% credible interval for the probability that the next flip will come up heads.

I believe:

I installed SciPy, and wrote:

from scipy.special import betainc, beta

numer = betainc(5257, 579112, 0.05)
denom = beta(5257, 579112)

print numer
print denom
print numer / denom

I get

./foo.py:6: RuntimeWarning: divide by zero encountered in double_scalars
  print numer / denom

That is not what I want. I believe this web page calculator says I want 0.008793824176376785502, which is plausible (close to 5257/(5257+579112)=0.00899602819451).

So: are my beliefs correct? How do I calculate this lower edge in Python?


Friend of mine answered this: I want the inverse cumulative distribution, i.e. the place on the X axis where the cumulative distribution is at 0.05.

In SciPy, the easiest thing is btdtri.

print btdtri(5257, 579112, 0.05)


| cite | improve this answer | |

Your Answer

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

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