Блог

Я теоретически вычисляю pmf на Python. вот код.

 >>> 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}.
в последней строке кода вычисляется вероятность для каждого из возможных результатов соответственно.

Я должен поместить 10 ‘a_coin’ в np.meshgrid (a_coin, …, a_coin), если я хочу вычислить pmf для подбрасывания 10 честных монет, что кажется скучным и неэффективным.

вопрос в том, есть ли более эффективный способ сделать это в 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)]))