Je calcule le pmf théoriquement en Python. voici le code.
>>> 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])
ce morceau de code simule 1 lancer de 3 pièces justes. les résultats possibles sont {0,1,2,3}. la dernière ligne de code calcule respectivement la probabilité pour chacun des résultats possibles.
Je dois mettre 10 'a_coin' dans np.meshgrid (a_coin, ..., a_coin) si je veux calculer le pmf pour lancer 10 pièces justes, ce qui semble ennuyeux et inefficace.
la question est, y a-t-il un moyen plus efficace de faire cela en python ou R?
3 Réponses :
Voici comment le faire dans R:
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
La fonction choose vous donne les coefficients binomiaux. Pour les transformer en probabilités, divisez simplement par leurs sommations:
> sapply(0:3, choose, n=3)/sum(sapply(0:3, choose, n=3)) [1] 0.125 0.375 0.375 0.125
Il ne semble pas que vous vouliez vraiment énumérer autant que calculer. Si vous avez besoin d'énumérer les résultats de 10 tirages binomiaux "équitables" successifs, vous pouvez utiliser combn 11 fois.
En utilisant les bibliothèques standard Python, vous pouvez obtenir des probabilités sous forme de nombres rationnels (c'est la solution exacte), par exemple
list(map(float, [Fraction(factorial(n), factorial(n - j)) * Fraction(1, factorial(j) * 2 ** n) for j in range(0, n + 1)]))
Cela pourrait être facilement converti en flottants, par exemple
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)]
Voici une solution fft basée sur numpy :
>>> 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])
Par exemple:
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]
Utilisez des formules mathématiques? Coeffs binomiaux etc ...