1. 程式人生 > 其它 >python實現二項分佈

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)