Can the scipy.stats.binom library return a value of "N" for a particular "k" and "p"

Question

I have what I hope is a simple question. I am trying to learn my way round SciPy as a hobby/spare-time activity (yes I really should get out more) and at the moment am playing with the scipy.stats.binom Binomial Probability functions (not helped by the fact it is 20+ years since I studied statistics).

So I have set myself the problem below.

"""
Binomial probability problem

The planet Klom is famous for two things. The superabundance of
dilithium crystals which are so common they are simply lying all
over the ground; and their notoriously poor quality as only 10%
of them are any good. Good and bad crystals are indistinguishable
to the naked eye.

Captain Kirk requires four good crystals to power his spaceship
and beams down to the surface. He can only carry 12 crystals back
to his ship.

What is the probability of him returning with 0, 1, 2, .. 12
good crystals.

Plot the probability

How many crystals should he bring back if he wants a better than 50,
90, 95 and 99% probability of having at least four working ones.
"""

import scipy.stats
import numpy
import matplotlib.pyplot

[N, p] = [12, 0.1]
rv = scipy.stats.binom(N, p)
x = numpy.arange(0, N+1)
pmf = rv.pmf(x)
for k in range(len(pmf)):
    print("{:3d} {:10.6f}".format(k, pmf[k]))

fig = matplotlib.pyplot.figure()
ax = fig.add_subplot(111)
ax.plot(x, pmf, 'bo')
ax.vlines(x, 0, pmf, lw=2)
ax.set_xlabel('Number of Good crystals')
ax.set_ylabel('Binomial PDF')
matplotlib.pyplot.show()

Now obviously I can solve the first part of the problem. But calling scipy.stats.binom(N, p) creates a frozen PMF in which the values of N and p are frozen and the probability computed for a given k.

How (other than by using a brute-force loop) can I work out the value of N which gives Kirk the 50%, 90%, etc. chance? (i.e. effectively working out N for a given k and p). I'm sure there must be some function that does it but I can't work out which.

Cheers.

Tim.

Answer 1

Since it is a discrete distribution we are dealing with here, we will need to use scipy.optimize

The first part is easy:

In [131]:

import scipy.stats as ss
import scipy.optimize as so
In [132]:

ss.binom.pmf(range(12), n=12, p=0.1)
Out[132]:
array([  2.82429536e-01,   3.76572715e-01,   2.30127770e-01,
         8.52325076e-02,   2.13081269e-02,   3.78811145e-03,
         4.91051484e-04,   4.67668080e-05,   3.24769500e-06,
         1.60380000e-07,   5.34600000e-09,   1.08000000e-10])

For the second part, we need to solve N form the equation PMF(4, N, 0.1)=0.5/0.1/0.05/0.01, by scipy.optimize, and we need the Nelder-Mead optimizer, which do not require information from the derivatives and thus suitable for the discrete problem we are facing.

In [133]:

for P in [0.5, 0.1, 0.05, 0.01]:
    N=int(so.fmin(lambda n, p: (ss.binom.pmf(4, int(n), 0.1)-p)**2, 40, args=(P,), disp=0)) 
    #the initial guess is 40, because 40*0.1=4
    print N,
    print ss.binom.pmf(4, N, p=0.1)
39 0.205887043441
67 0.100410451946
81 0.0498607360095
107 0.00999262321 #need to pick 107 crystals in order to be sure 99% of time that there are at least 4 usable crystals.