scipy.stats.binom库可以返回特定"k“和"p”的"N“值吗？

Question

我所希望的是一个简单的问题。我正在努力学习如何绕过SciPy，以此作为业余爱好/业余活动(是的，我真的应该多出去)，目前我正在玩scipy.stats.binom二项式概率函数(这对我学习统计学已经是20+年的事实没有帮助了)。

所以我给自己设置了下面的问题。

"""
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()

显然，我可以解决问题的第一部分。但是调用scipy.stats.binom(N, p)会创建一个冻结的PMF，其中N和p的值被冻结，并计算给定k的概率。

除了使用蛮力循环之外，我如何计算给柯克50%、90%等机会的N值呢？(即有效计算给定k和p的N)。我肯定一定有什么功能可以做，但我想不出是哪一个。

干杯。

蒂姆。

Answer 1

由于我们在这里处理的是一个离散的分布，所以我们需要使用scipy.optimize。

第一部分很简单：

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])

对于第二部分，我们需要用scipy.optimize求解方程PMF(4，N，0.1)=0.5/0.1/0.0 5/0.0 1，并且需要Nelder-Mead优化器，它不需要导数的信息，适合于我们所面临的离散问题。

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.