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 could use a beta distribution with alpha=5257 and beta=579112.
- I want the value of the cumulative distribution at 0.05.
- I could calculate that in Python with SciPy using the betainc and beta special functions.
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
1.0 0.0 ./foo.py:6: RuntimeWarning: divide by zero encountered in double_scalars print numer / denom inf
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?