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)