Prawdopodobieństwo, że wartość zmiennej losowej \(X\) znajdzie się w zadanym przedziale, można wyznaczyć również z symulacji. Możliwość ta przydaje się zarówno do zmiennych ciągłych jak i dyskretnych, gdy nie znamy analitycznej postaci dystrybuanty lub chcemy po prostu zweryfikować wyniki otrzymane analitycznie.
Przykładowo, zbadajmy ile wynosi \(P(X<0)\) dla \(X \sim N(\mu=20, \sigma=10)\) czyli prawdopodobieństwo, że wartość wylosowana z rozkładu normalnego o \(\mu=20\) i \(\sigma=10\) będzie ujemna.
import numpy, scipy.stats
ile_wszystkich = 100000 # tyle liczb wygenerujemy
rozklad = scipy.stats.norm(20,10)
# generujemy N liczb z rozkładu normalnego o zadanych parametrach
tablica = rozklad.rvs(size=ile_wszystkich)
# zliczamy ile jest liczb ujemnych
ile_ujemnych = 0
for x in tablica :
if x < 0 :
ile_ujemnych += 1
# lub nieco krócej
ile_ujemnych = numpy.sum(tablica < 0)
# prawdopodobieństwo (uwaga na dzielenie liczb całkowitych)
p_symulacja = float(ile_ujemnych) / float(ile_wszystkich)
# możemy to porównać z prawdopodobieństwem analitycznym
p_analityczne = rozklad.cdf(0) # dla rozkładów ciągłych P(X<0) równa się P(X≤0)
Analogicznie, zbadajmy ile wynosi \(P(Y \ge 2)\) dla \(Y \sim B(n=10, p=0.1)\) (rozkład dwumianowy) czyli prawdopodobieństwo, że otrzymamy co najmniej dwa sukcesy na dziesięć prób, jeśli prawdopodobieństwo pojedynczego sukcesu to \(0.1\).
import numpy, scipy.stats
ile_wszystkich = 100000 # tyle liczb wygenerujemy
rozklad = scipy.stats.binom(10,0.1)
# generujemy N liczb z rozkładu dwumianowego o zadanych parametrach
tablica = rozklad.rvs(size=ile_wszystkich)
# zliczamy ile jest liczb większych lub równych 2
ile_wr2 = 0
for x in tablica :
if x >= 2 :
ile_wr2 += 1
# lub nieco krócej
ile_wr2 = numpy.sum(tablica >= 2)
# prawdopodobieństwo (uwaga na dzielenie liczb całkowitych)
p_symulacja = float(ile_wr2) / float(ile_wszystkich)
# możemy to porównać z prawdopodobieństwem analitycznym
# równym P(Y>=2) czyli P(Y>1) czyli 1 - P(Y<=1)
p_analityczne = 1 - rozklad.cdf(1)
# można też policzyć P(Y<=1) = P(Y=0) + P(Y=1)
p_analityczne = 1 - (rozklad.pmf(0) + rozklad.pmf(1))
autor: Piotr Różański, ostatnia modyfikacja: 10.04.2016