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 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
I get
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?
Best Answer
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.
gives