In [29]:
import random
import matplotlib.pyplot as plt
from scipy import stats
from math import ceil
import numpy as np
Geometrické rozdělení na koleně¶
Budeme generovat náhodnou veličinu s geometrickým rozdělením přímo podle definice. A když jich tak vygenerujeme hodně, odhadneme pravděpodobnost jevů $X = 1$, $X=2$, atd. -- prostě tak, že spočítáme, v jakém procentu případů jsme takové číslo vygenerovali.
In [34]:
def plot_geom(p = 1/6, n_samples = 100, max_k = 20):
"""
p je pravděpodobnost úspěchu
n_samples je počet vzorků (kolikrát čekáme na šestku)
max_k je max. hodnota zobrazená v grafu
"""
# Jedna hodnota geometrického rozdělení
def geo_sample(p):
count = 1
while random.random() > p: # skončíme s pravděpodobností p
count += 1
return count
samples = [geo_sample(p) for _ in range(n_samples)] # naházíme kostkami a zapíšeme si dobu čekání na šestku
k_vals = list(range(1, max_k+1)) # pro tyhle hodnoty budeme kreslit graf
freq = [ samples.count(k)/n_samples for k in k_vals ] # spočteme počty výskytů a vydělíme počtem opakování -- odhad pravděpodobnosti
prob = [(1-p)**(k-1) * p for k in k_vals ] # teoretická pravděpodobnost
# Plot results
plt.figure(figsize=(8, 5))
plt.bar(k_vals, freq, width=0.7, label="Naměřené", alpha=0.6)
plt.plot(k_vals, prob, 'ro', markersize=4, label="Teoretické")
plt.xlabel("k (Počet pokusů do úspěchu)")
plt.ylabel("Pravděpodobnost")
plt.xticks(k_vals)
plt.title(f"Geo(p={p:.3}) ")
plt.legend()
plt.grid(True)
plt.show()
In [31]:
random.seed(10)
plot_geom()
Ilustrace zacházení s náhodnými veličinami ve scipy¶
In [36]:
X = stats.geom(p=1/6) # X je náhodná veličina -- generátor náhodných čísel s danými pravděpodobnostmi jednotlivých výstupů
In [42]:
xvals = X.rvs(100) # vygenerujeme 100 hodnot
list(xvals).count(1)
Out[42]:
22
In [38]:
X.pmf(1) # necháme si spočítat pravděpodobnost vygenerování 1
Out[38]:
0.16666666666666666
In [45]:
X.mean() # ideální střední hodnota
Out[45]:
6.0
In [46]:
xvals.mean() # naměřená střední hodnota -- průměr hodnot z xvals (využíváme tady efektivitu knihovny numpy)
Out[46]:
5.66
In [35]:
def plot_any(X, text = "X", n_samples = 100, min_k = 1, max_k = 20, scale=20, filename = None, save_format = None):
"""
X je scipy náhodná veličina
n_samples je počet vzorků (kolikrát čekáme na šestku)
max_k je max. hodnota zobrazená v grafu
"""
samples = list(X.rvs(size=n_samples)) # "Naházíme" několik hodnot náhodné veličiny
k_vals = list(range(min_k, max_k+1)) # Spočteme počty výskytů
freq = [ samples.count(k)/n_samples for k in k_vals ]
prob = [ X.pmf(k) for k in k_vals ]
p_max = ceil(max(max(prob),max(freq))*scale)
p_vals = [k/scale for k in range(p_max+1)]
# Plot results
plt.figure(figsize=(8, 5))
plt.bar(k_vals, freq, width=0.7, label="Naměřená", alpha=0.6)
plt.plot(k_vals, prob, 'ro', markersize=4, label="Teoretická")
plt.xlabel("k")
plt.ylabel("Pravděpodobnost")
plt.xticks(k_vals)
plt.yticks(p_vals)
plt.title(text + " -- teoretická vs. empirická (naměřená) pravděpodobnost")
plt.legend()
plt.grid(True)
if save_format in ['png', 'jpg', 'pdf']:
if filename is None:
filename = text
full_filename = f"{filename}.{save_format}"
plt.savefig(full_filename, dpi=300, bbox_inches='tight')
print(f"Plot saved as {full_filename}")
plt.show()
Generování obrázků pro různá rozdělení pomocí scipy¶
In [33]:
np.random.seed(10) # scipy.stats používá generátor z numpy, ne ten balíčku random ve standardním Pythonu
plot_any(stats.geom(p=1/6), text="Geo(1/6)", n_samples = 100, scale=100, min_k = 1, max_k = 20, filename="Geo(0.167)", save_format = "png")
Plot saved as Geo(0.167).png
In [ ]:
In [6]:
np.random.seed(37)
plot_any(stats.binom(n=10, p=1/6), text="Bin(10,1/6)", n_samples = 20, min_k = 0, max_k = 10, filename="Bin(10,0.167)", save_format = "png")
Plot saved as Bin(10,0.167).png
In [7]:
np.random.seed(43)
plot_any(stats.poisson(10/6), text="Poi(10/6)", n_samples = 100, min_k = 0, max_k = 10, filename="Poi(1.67)", save_format = "png")
Plot saved as Poi(1.67).png
Poissonova aproximace¶
Zafixujeme $\lambda > 0$ a $k$. Pak pravděpodobnost, že náhodná veličina s rozdělením $Bin(n,\lambda/n)$ nabude hodnotu $k$ se limitně blíží pravděpodobosti, že se $k$ rovná Poissonova náhodná veličina s parametrem $\lambda$.
Proč je to zajímavé: pro málo častý jev (v daném časovém intervalu dostaneme email/SMS/dojde k nehodě/...) platí, že když dvakrát zmenšíme časový úsek, tak se dvakrát zmenší pravděpodobnost dané události. A zároveň dvakrát zvýší počet "pokusů". Přesněji, pokud si časový úsek rozdělíme na $n$ dílů, tak máme $n$ pokusů a $n$-krát nižší pravděpodobnost, což vede k binomickému rozdělení s parametry $Bin(n, \lambda/n)$.
In [115]:
def cmp_bin_Pois(n=10, p = 1/6, n_samples = 100, min_k = 0, max_k = 15, a = None, b = None, save_format=None, filename=None):
"""
n_samples je počet vzorků (kolikrát čekáme na šestku)
max_k je max. hodnota zobrazená v grafu
"""
X = stats.binom(n,p)
Y = stats.poisson(n*p)
k_vals = list(range(min_k, max_k+1))
p_X = [ X.pmf(k) for k in k_vals ]
p_Y = [ Y.pmf(k) for k in k_vals ]
# Plot results
plt.figure(figsize=(10, 5))
plt.plot(k_vals, p_X, 'ro', markersize=4, label=a)
plt.plot(k_vals, p_Y, 'bx', markersize=4, label=b)
plt.xlabel("k")
plt.ylabel("Pravděpodobnost")
plt.xticks(k_vals)
plt.title(f"Poissonova approximace: {a} vs. {b}",)
plt.legend()
plt.grid(True)
if save_format in ['png', 'jpg', 'pdf']:
if filename is None:
filename = f"bin_pois_approx_10_1:6"
full_filename = f"{filename}.{save_format}"
plt.savefig(full_filename, dpi=300, bbox_inches='tight')
print(f"Plot saved as {full_filename}")
plt.show()
In [212]:
cmp_bin_Pois(a = "Bin(10,1/6)", b = "Poi(10/6)", save_format="png")
Plot saved as bin_pois_approx_10_1:6.png
In [164]:
cmp_bin_Pois(n=100, p=1/60, a = "Bin(100,1/60)", b = "Poi(100/60)", filename="bin_pois_approx_100_1:60", save_format="png")
Plot saved as bin_pois_approx_100_1:60.png
In [ ]: