Python中的二项式测试,用于非常大的数字

Question

我需要在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兼容代码.

非常感谢!

Answer 1

编辑添加此评论:请注意,正如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

BTW:适用于scipy 0.7.2,以及当前的0.8 dev.

Answer 2

任何看起来像的解决方案comb(n, k) * 0.5**k * 0.5**(n-k)都不适用于大型解决方案n.在大多数(所有？)平台上,Python float可以存储的最小值大约为2** - 1022.对于大n-k或大k,右侧将四舍五入为0.同样,梳子(n,k)可以变得如此之大,以至于它不适合浮子.

更稳健的方法是将概率密度函数计算为累积分布函数中两个连续点之间的差异,这可以使用正则化的不完全β函数计算(参见SciPy的"特殊函数"包).数学:

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

我没有测试过代码,但这应该给你一般的想法.