python实现二项分布

二项分布的定义

在概率论和统计学中，二项分布是n个独立的成功/失败试验中成功的次数的离散概率分布，其中每次试验的成功概率为p

如果随机变量X服从参数为n和p的二项分布，我们记为X~B(n,p)

n次试验中正好得到k次成功的概率由概率质量函数给出：\(P\{X=k\}={C_n^k}{p^k}{(1-p)^{n-k}}\), 其中 \({C_n^k}=\frac{n!}{k!(n-k)!}\)

代码实现

from functools import reduce

def factorial(n):
    """计算n的阶乘,即n!"""
    if n == 0:
        return 1
    return reduce(lambda x,y: x*y, range(1, n+1))

def Cnk(n, k):
    """计算从n个样本中取出k个样本，共有多少种组合"""
    return int(factorial(n)/(factorial(k) * factorial(n-k)))

def binomial(n, p, k):
    """当X服从二项分布B(n,p)时，X取值为k的概率"""
    return Cnk(n, k) * (p**k) * ((1-p)**(n-k))

然而在实际中，如果n取比较大的值，那么上述代码在运行时就会产生数值溢出的问题，为了解决这个问题，将目标改成求概率质量函数的对数值，即 \(log(P\{X=k\})\)

代码

import numpy as np

def log2_sum(a, b):
    """计算log2(a)+log2(a+1)+...+log2(b)"""
    return sum([np.log2(i) for i in range(a, b+1)])

def log2_binomial(n, p, k):
    log = np.log2
    return log2_sum(k+1, n) - log2_sum(1, n-k) + k * log(p) + (n-k) * log(1-p)