Jak wyznaczyć prawdopodobieństwo otrzymania zadanych wartości zmiennej losowej?

dla dowolnej zmiennej losowej, korzystając z symulacji

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))