Popis problému

lambda_ = 10

shape_ = 2 scale_ = 0.25

Nastavení a importy

import numpy as np from scipy.stats import expon, poisson, gamma, chi2, chisquare from matplotlib import pyplot as plt

Úlohy

Vytváříme trajektorii Poissonova procesu s obsluhou, kde intervaly mezi příchody jsou exponenciálně rozdělené s parametrem míry lambda_ a délky obsluhy jsou gama rozdělené s parametry tvaru shape_ a měřítka scale_.

def create_trajectory(T: float, lambda_, shape_, scale_): proces = [[0, 0]] t = 0 in_scale=1/lambda_ while t < T: t_delta = expon.rvs(scale=in_scale) if t + t_delta > T: break proces.append([t+t_delta, 1]) t_obsluhy = gamma.rvs(shape_, scale=scale_) if t+t_delta+t_obsluhy <= T: proces.append([t+t_delta+t_obsluhy, -1]) t += t_delta proces = np.array(proces) return proces[np.argsort(proces[:, 0])]

proces_1 = create_trajectory(10, lambda_, shape_, scale_)

Průběh trajektorie graficky znázorněme: na ose x znázorňujeme hodnoty z prvního sloupce pole proces_1, což jsou časové okamžiky, kdy se vyskytly jednotlivé události (příchody a konce obsluh). Na ose y znázorňujeme kumulativní součet hodnot z druhého sloupce pole proces_1, které představují změny v Poissonově procesu v čase (počet příchodu minus počet konců obsluh), pomocí funkce z knihovny numpy - cumsum https://numpy.org/doc/stable/reference/generated/numpy.cumsum.html

plt.plot(proces_1[:, 0], np.cumsum(proces_1[:, 1])) plt.title("Nahodný proces M∣G∣∞: 10 sec") plt.xlabel('Čas') plt.ylabel('Počet zákazníků')

Nastavíme požadovane parametry:

N = 500 T = 100

Vytváříme N náhodných trajektorií pomocí funkce create_trajectory s parametry (T, lambda_, shape_, scale_).

processes_results = np.zeros(N) for i in np.arange(N): trajectory = create_trajectory(T, lambda_, shape_, scale_) processes_results[i] = np.sum(trajectory[:, 1])

trajectory

Každá trajektorie je reprezentována jako dvousloupcová matice, kde první sloupec jsou časové okamžiky, kdy dochází ke změnám v počtu zákazníků, druhý sloupec obsahuje "1" a "-1" (změny v počtu klientů)

Graficky znázorněme simulované trajektorie.

plt.plot(np.arange(N), processes_results) plt.title("Simulace 500 nezávislých trajektorií M∣G∣∞: 100 sek ") plt.xlabel('Pokus') plt.ylabel(f'Počet zákazníků po {T} sekund')

Trajektorie seskupíme podle počtu zákazníků ve frontě po 100 sekundách.

Ni = np.arange(int(np.max(processes_results)+1)) for i in Ni: Ni[i] = np.sum(processes_results == i) print(sum(Ni), Ni) Ni_distribution = Ni / np.sum(Ni) plt.bar(np.arange(len(Ni_distribution)), Ni_distribution) plt.title("Rozdělení počtů zákazníků za 100 sekund") plt.xticks(np.arange(len(Ni_distribution))) plt.xlabel('Číslo zákazníků') plt.ylabel('Pravděpodobnost')

Vidíme, že tato distribuce je velmi podobná normálnímu rozdělení.

poisson_intenzita = lambda_ * shape_ * scale_ poisson.mean(poisson_intenzita)

Vypočítáme a graficke znázorněme distribuce Poissonova rozdělení s určenou intenzitou (5).

poisson_probs = np.zeros(int(np.max(processes_results)+1), dtype=float) for i in np.arange(len(Ni_distribution)): poisson_probs[i] = poisson.pmf(i, poisson_intenzita) poisson_probs[-1] = 1 - np.sum(poisson_probs[:-1]) print(poisson_probs) plt.bar(np.arange(len(poisson_probs)), poisson_probs) plt.title("Rozdělení Poisson s intenzitou 5") plt.xticks(np.arange(len(Ni_distribution))) plt.xlabel('Číslo zákazníků') plt.ylabel('Pravděpodobnost')

Testování

Ni_casted = Ni.copy() p_casted = poisson_probs * N X_casted = np.arange(int(np.max(processes_results)+1)) p_casted

while np.min(p_casted) < 5: order = np.argsort(p_casted) p_casted[order[1]] += p_casted[order[0]] Ni_casted[order[1]] += Ni_casted[order[0]] p_casted = p_casted[~np.isin(np.arange(p_casted.shape[0]), order[0])] Ni_casted = Ni_casted[~np.isin(np.arange(Ni_casted.shape[0]), order[0])] X_casted = X_casted[~np.isin(np.arange(X_casted.shape[0]), order[0])] print(np.concatenate((p_casted.reshape(-1, 1), Ni_casted.reshape(-1, 1)), axis=1)) plt.bar(X_casted, p_casted, alpha=0.5, label = 'Poisson') plt.bar(X_casted, Ni_casted, alpha=0.5, label = 'Proces') plt.xticks(X_casted) plt.title("Porovnání limitního rozdělení (Poisson) s naším procesem (Simulace)") plt.xlabel('Číslo zakazniků') plt.ylabel('Hodnoty') plt.legend()

Spočítáme testovou statistiku.

test_stat = 0 for i, p in enumerate(p_casted): test_stat += (Ni_casted[i] - p)**2 / (p) print("Testová statistika:" ,test_stat)

Dále, abychom otestovali hypotézu, musíme najít kritickou hodnotu. K tomu použijeme knihovnu scipy.stats k výpočtu kritické hodnoty pro test chí-kvadrát na hladině významnosti alfa = 0,05. Kritická hodnota se vypočítá pomocí funkce stats.chi2.ppf(). https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html Tato funkce se používá k určení kritické hodnoty pro test chí-kvadrát na dané hladině významnosti (5%) a odpovídajících stupňů volnosti.

print("Počet stupňů volnosti:" ,p_casted.shape[0] - 1)

# kritický obor alpha = 0.05 krit_hodnota = chi2.ppf(1-alpha, p_casted.shape[0] - 1) print("Kritická hodnota:" ,krit_hodnota)

test_result = test_stat > krit_hodnota print(f"test_stat = {np.round(test_stat, 3)}, kritická_hodnota = {np.round(krit_hodnota, 3)}, test_stat {'>' if test_result else '<'} krit_hodnota =>") print(f"=> hypotézu, že skutečné hodnoty pravděpodobností jsou Poisson({poisson_intenzita}), {'ZAMITNEME' if test_result else 'SCHVÁLÍME'}")

Pokud Testová statistika < Kritická hodnota =>

Určení p-hodnoty testu

K určení p-hodnoty také použijeme knihovnu scipy.stats.

Najdeme p-hodnotu pomocí scipy.stats.chi2.sf - survival function (1 - cdf)

print(f"p-hodnota = {chi2.sf(test_stat, p_casted.shape[0] - 1)}")

p-hodnota je větší než hladina významnosti. Platí tedy, že čím výšší p-hodnota testu je, tím větší nám tento test indikuje pravděpodobnost, že platí nulová hypotéza.

Pro jistotu můžeme naše výsledky porovnat s výsledky získanými pomocí knihovní funkce - scipy.stats.chisquare(), která vypočítává jednostranný chí-kvadrát test.

# pro jistotu chisquare(Ni_casted, p_casted)

.css-15w88e5{color:var(--chakra-colors-fg-neutral-primary);font-weight:inherit;letter-spacing:-0.09px;}Popis problému

Nastavení a importy

Úlohy

Testování

Určení p-hodnoty testu

Popis problému