我需要在Python中进行二项式检验，可以计算10000阶数的'n'个数字。
我已经使用scipy.misc.comb实现了一个快速的binomial_test函数，然而，它在n = 1000左右就受到了很大的限制，我猜测是因为在计算阶乘或组合数本身时达到了最大可表示的数字。这是我的函数：
from scipy.misc import comb
def binomial_test(n, k):
"""Calculate binomial probability
"""
p = comb(n, k) * 0.5**k * 0.5**(n-k)
return p
我该如何使用本地的Python（或numpy、scipy等）函数来计算二项式概率？如果可能的话，我需要与scipy 0.7.2兼容的代码。
非常感谢！
编辑添加此评论：请注意，正如Daniel Stutzbach所提到的，"二项式检验"可能不是原帖作者想要的（尽管他使用了这个表达）。他似乎是在问二项分布的概率密度函数，而这不是我下面建议的内容。
你尝试过scipy.stats.binom_test吗？
rbp@apfelstrudel ~$ python
Python 2.6.2 (r262:71600, Apr 16 2009, 09:17:39)
[GCC 4.0.1 (Apple Computer, Inc. build 5250)] on darwin
Type "help", "copyright", "credits" or "license" for more information.
>>> from scipy import stats
>>> print stats.binom_test.__doc__
Perform a test that the probability of success is p.
This is an exact, two-sided test of the null hypothesis
that the probability of success in a Bernoulli experiment
is `p`.
Parameters
----------
x : integer or array_like
the number of successes, or if x has length 2, it is the
number of successes and the number of failures.
n : integer
the number of trials. This is ignored if x gives both the
number of successes and failures
p : float, optional
The hypothesized probability of success. 0 <= p <= 1. The
default value is p = 0.5
Returns
-------
p-value : float
The p-value of the hypothesis test
References
----------
.. [1] http://en.wikipedia.org/wiki/Binomial_test
>>> stats.binom_test(500, 10000)
4.9406564584124654e-324
添加文档链接的小修改: http://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.binom_test.html#scipy.stats.binom_test
顺便说一下：在scipy 0.7.2上可以使用，也可以在当前的0.8 dev版本上使用。
comb(n, k) * 0.5**k * 0.5**(n-k) 的解决方案都不适用于大的 n。在大多数（所有？）平台上，Python浮点数可以存储的最小值约为2 ** -1022。对于大的n-k或k，右侧将被舍入为0。同样，comb(n, k)可以增长到无法适应浮点数。更健壮的方法是计算概率密度函数，作为累积分布函数中两个连续点之间的差异，可以使用正则化不完全beta函数进行计算（请查看SciPy的“special functions”包）。数学上：pdf(p, n, k) = cdf(p, n, k) - cdf(p, n, k-1)
另一种选择是使用正态分布近似，对于大的n来说非常准确。如果速度是一个问题，那么这可能是最好的选择：
from math import *
def normal_pdf(x, m, v):
return 1.0/sqrt(2*pi*v) * exp(-(x-m)**2/(2*v))
def binomial_pdf(p, n, k):
if n < 100:
return comb(n, k) * p**k * p**(n-k) # Fall back to your current method
return normal_pdf(k, n*p, n*p*(1.0-p))
我还没有测试过这段代码，但它应该能让你大致明白。
comb(n, k) * p**k * p**(n-k) 应该是 comb(n, k) * p**k * (1-p)**(n-k)。 - minillinimGMPY 还支持扩展精度的浮点数计算。例如：
>>> from gmpy import *
>>>
>>> def f(n,k,p,prec=256):
... return mpf(comb(n,k),prec) * mpf(p,prec)**k * mpf(1-p,prec)**(n-k)
...
>>> print(f(1000,500,0.5))
0.0252250181783608019068416887621024545529410193921696384762532089115753731615931
>>>
我会研究GNU多精度库(gmpy)，它允许您执行任意精度的计算：您可能可以这样做：
comb(n, k, exact=1)/2**k/2**(n-k)
但使用gmpy的长整数。
实际上，如果您使用精确的整数计算，您可以轻松地达到n=10000 用于组合部分; 为此，您必须使用：
comb(n, k, exact=1)
而不是浮点近似值 comb(n, k)，会导致溢出。
但正如帖子作者所指出的那样，返回的（长）整数可能太长无法乘以浮点数！
此外，我们很快遇到另一个问题：0.5**1000=9.3…e-302 已非常接近浮点下溢…
总之：如果您真的需要所有 k 对于 n~10,000 的精确结果，则需要使用与原始帖子中的公式不同的方法，该公式受双精度浮点算术的限制。 使用上面提到的 gmpy 可能是解决方案（未经测试！）。
并非特定针对 Python 的解决方案，但如果您可以处理小的分数误差，则可以尝试使用 Stirling 公式逼近 n!：
comb(n, k) = n!/(k! * (n-k)!), 这里 n! 在 n 很大时近似于 sqrt(2*Pin)(n/e)^n。
对于 n＞1000，分数误差应该非常小。
对于大 n 的概率计算，使用对数进行中间结果的计算：
log p = log(comb(n, k)) - n * log(2)
p = exp(log(p))
# This imports the array function form numpy
from numpy import array
# the following defines the factorial function to be used in the binomial commands/
# n+1 is used in the range to include the nth term
def factorial (n):
f=1
for x in range(1,n+1):
f=f*(x)
return f
# The follwong calculates the binomial coefficients for given values of n & k
def binomial (n,k):
b=1
b=(factorial(n)/(factorial(k)*factorial(n-k)))
return int(b)
# the following lines define the pascal triangle , and print it out for 20 rows./
# in order to include nth term, the n +1 term needs to be in the range. The commands/
# append the next binomial coeficiant to a raw first and then append rows to the triangle/
# and prints a 20 row size pascal triangle
def pascal(T):
triangle=[]
for n in range(T):
r=[]
for k in range(n+1):
r.append(binomial(n,k))
triangle.append(r)
return triangle
for r in pascal(20):
print((r))