我在Python中理论上计算了pmf。这是代码。
>>> a_coin = np.array([0,1])
>>> three_coins = np.array(np.meshgrid(a_coin,a_coin,a_coin)).T.reshape(-1,3)
>>> heads = np.sum(three_coins, axis = 1)
>>> df = pd.DataFrame({'heads': heads, 'prob': 1/8})
>>> np.array(df.groupby('heads').sum()['prob'])
array([0.125, 0.375, 0.375, 0.125])
这段代码模拟了1折3公平硬币。可能的结果是{0,1,2,3}。最后一行代码分别计算每个可能结果的概率。
我必须在np.meshgrid(a_coin,...,a_coin)中放入10'a_coin',如果我想计算投掷10个公平硬币的pmf,这看起来很无聊且效率低下。
问题是,在python或R中有更有效的方法吗?
以下是如何在R中执行此操作:
> sapply(0:3, choose, n=3)/sum(sapply(0:3, choose, n=3))
[1] 0.125 0.375 0.375 0.125
choose
函数为您提供二项式系数。将它们变成概率只需除以它们的总和:
sapply(0:10, choose, n=10)
[1] 1 10 45 120 210 252 210 120 45 10 1
sapply(0:10, choose, n=10)/ sum( sapply(0:10, choose, n=10))
[1] 0.0009765625 0.0097656250 0.0439453125 0.1171875000 0.2050781250 0.2460937500 0.2050781250
[8] 0.1171875000 0.0439453125 0.0097656250 0.0009765625
看起来你真的不想计算这么多。如果你需要枚举10次连续“公平”二项式抽签的结果,那么你可以使用combn
11次。
这是一个基于fft
的numpy
解决方案:
import numpy as np
from scipy import fftpack
def toss(n=10, p=0.5):
t1 = np.zeros(fftpack.next_fast_len(n+1))
t1[:2] = 1-p, p
f1 = fftpack.rfft(t1)
c1 = f1[1:(len(t1) - 1) // 2 * 2 + 1].view(f'c{2*t1.itemsize}')
c1 **= n
f1[::(len(t1) + 1) // 2 * 2 - 1] **= n
return fftpack.irfft(f1)[:n+1]
例如:
>>> toss(3)
array([0.125, 0.375, 0.375, 0.125])
>>> toss(10)
array([0.00097656, 0.00976562, 0.04394531, 0.1171875 , 0.20507813,
0.24609375, 0.20507813, 0.1171875 , 0.04394531, 0.00976562,
0.00097656])
使用Python标准库,您可以将概率作为有理数(这是精确解),例如
from fractions import Fraction
from math import factorial
n=30
[Fraction(factorial(n), factorial(n - j)) * Fraction(1, factorial(j) * 2 ** n) for j in range(0, n + 1)]
这可以很容易地转换成浮点数,例如
list(map(float, [Fraction(factorial(n), factorial(n - j)) * Fraction(1, factorial(j) * 2 ** n) for j in range(0, n + 1)]))