NI-VSM-DU-3

Data a parametry

Reprezentant - Bogdan Buliakov

Parametry úlohy byly spočtěny následovně:

K = den narození reprezentanta skupiny - 23 L = počet písmen v příjmení reprezentanta - 8 (BULIAKOV)

Podle uvedeného vzorce jsme definovali datový soubor č. 1, jako soubor 013.txt a datový soubor č. 2, jako 014.txt

K = 23 L = 8 X = ((K*L*23) % (20)) + 1 Y = ((X + ((K*5 + L*7) % (19))) % (20)) + 1 print('X =', X) print('Y =', Y)

X = 13
Y = 14

Nastavení

import scipy import numpy as np from scipy import stats from matplotlib import pyplot as plt %matplotlib inline

Funkce read_text_from_file se používá ke čtení textu ze souboru

Funkce abs_freq se používá ke výpočtu absolutní četnosti jednotlivých znaků

Funkce draw_abs_freq se používá ke grafickému znázornění absolutní četnosti jednotlivých znaků

class VSM: def __init__(self): super().__init__() self.text = "" self.text_code = [] self.letters_codes = [] self.letters = [] def read_text_from_file(self, file_str: str): with open(f"{file_str}", 'r') as f: self.text = f.readlines()[1] self.text_code = list([ord(x) for x in self.text]) def abs_freq(self): self.letters_codes, self.letters_counts = np.unique(self.text_code, return_counts=True) self.letters_abs_freq = self.letters_counts self.letters = [chr(int(c)) for c in self.letters_codes] self.alphabet_len = self.letters_codes.shape[0] return list(zip(self.letters, self.letters_abs_freq)) def draw_abs_freq(self): fig, ax = plt.subplots() ax.bar(np.arange(len(self.letters)), self.letters_abs_freq) fig.tight_layout() plt.xticks(ticks=np.arange(len(self.letters)), labels=self.letters) plt.title("Absolutní četnosti jednotlivých znaků (symbolů včetně mezery)") plt.xlabel("Znak") plt.ylabel("Absolutní četnost")

Z obou datových souborů načteme texty k analýze. Pro každý text zvlášť zjistime absolutní četnosti jednotlivých znaků (symbolů včetně mezery), které se v textech vyskytují pomocí výše popsaných funkcí.

Výsledky pro první text

print(f"Text 1: soubor {X}") vsm1 = VSM() vsm1.read_text_from_file(f"{str(X).zfill(3)}.txt") vsm1.abs_freq() vsm1.draw_abs_freq() print(f"Absolutní četnosti jednotlivých znaků (symbolů včetně mezery)") for item in vsm1.abs_freq(): print(f"'{item[0]}': {item[1]}")

Text 1: soubor 13
Absolutní četnosti jednotlivých znaků (symbolů včetně mezery)
' ': 1359
'a': 513
'b': 83
'c': 180
'd': 256
'e': 819
'f': 154
'g': 113
'h': 323
'i': 471
'j': 10
'k': 20
'l': 225
'm': 195
'n': 433
'o': 420
'p': 108
'q': 9
'r': 371
's': 387
't': 516
'u': 150
'v': 92
'w': 110
'x': 11
'y': 157
'z': 4

Výsledky pro druhý text

print(f"Text 2: soubor {Y}") vsm2 = VSM() vsm2.read_text_from_file(f"{str(Y).zfill(3)}.txt") vsm2.abs_freq() vsm2.draw_abs_freq() print(f"Absolutní četnosti jednotlivých znaků (symbolů včetně mezery)") for item in vsm2.abs_freq(): print(f"'{item[0]}': {item[1]}")

Text 2: soubor 14
Absolutní četnosti jednotlivých znaků (symbolů včetně mezery)
' ': 1139
'a': 424
'b': 78
'c': 121
'd': 236
'e': 652
'f': 134
'g': 98
'h': 337
'i': 328
'j': 7
'k': 35
'l': 170
'm': 113
'n': 356
'o': 386
'p': 76
'q': 4
'r': 335
's': 334
't': 476
'u': 133
'v': 50
'w': 125
'x': 4
'y': 92
'z': 2

1. Za předpokladu výše odhadněte matici přechodu markovského řetězce pro první text. Pro odhad matice přechodu vizte přednášku 17. Odhadnuté pravděpodobnosti přechodu vhodně graficky znázorněte, např. použitím heatmapy.

Funkce calc_transition_matrix se používá k odhadu matice přechodu markovského řetězce pro první text. Jestli předpokládejme, že první text je vygenerován z homogenního markovského řetězce s diskrétním časem, tak budeme odhadovat matice přechodu pomocí četností přechodů (Přednáška 17):

Poprve vytvoříme nulovou matici P_matice s rozměry [self.alphabet_len] x [self.alphabet_len] (self.alphabet_len - tato proměnná představuje počet unikátních písmen (znaků) v daném textu, které byly vypočtena ve funkci abs_freq)

Pak vytvoříme slovník letters_dict, který přiřazuje každému písmenu v abecedě číslo od 0 do [self.alphabet_len] - 1.

Potom je krok s procházením textu a aktualizováním hodnot v P_Matice na základě počtu přechodů z jednoho písmena na druhé.

Na konci je normalizace hodnot v P_matice tak, aby součet hodnot v každém řádku byl roven 1.

Funkce draw_transition_matrix se používá ke grafickému znázornění odhadnuté pravděpodobnosti přechodu, použitím heatmapy.

def calc_transition_matrix(vsm): P_matice = np.zeros((vsm.alphabet_len, vsm.alphabet_len)) letters_dict = {} for i, lt in enumerate(vsm.letters): letters_dict[lt] = i prev_ltr = vsm.text[0] for ltr in vsm.text[1:]: P_matice[letters_dict[prev_ltr], letters_dict[ltr]] += 1 prev_ltr = ltr P_matice = P_matice / np.sum(P_matice, axis=1).reshape(-1, 1) return P_matice def draw_transition_matrix(vsm, P_matice): fig, ax = plt.subplots(dpi=120) P_show = ax.imshow(P_matice) plt.xticks(ticks=np.arange(len(vsm.letters)), labels=vsm.letters) plt.yticks(ticks=np.arange(len(vsm.letters)), labels=vsm.letters) plt.title("Matice přechodu") ax.xaxis.tick_top() ax.xaxis.set_label_position('top') ax.set_xlabel("Cílový znak", loc='left') ax.set_ylabel("Zdrojový znak", loc='top') cbar = fig.colorbar(P_show)

Odhadneme matici přechodu markovského řetězce pro první text a odhadnuté pravděpodobnosti přechodu vhodně graficky znázorneme pomocí výše popsaných funkcí.

P_matice = calc_transition_matrix(vsm1) draw_transition_matrix(vsm1, P_matice)

Součet v každem řádku pro kontrolu:

np.sum(P_matice, axis=1)

2. Na základě matice z předchozího bodu najděte stacionární rozdělení 𝜋 tohoto řetězce pro první text.

Výsledky pro první text

Funkce stationary_distribution se používá ke výpočtu stacionárního rozdělení (Přednáška 14):

def stationary_distribution(vsm, P_matice): # NI-VSM-Lec-14-Handout.pdf stat_dist = scipy.linalg.null_space((P_matice - np.eye(vsm.alphabet_len)).T).T[0] stat_dist = stat_dist / np.sum(stat_dist) return stat_dist

pi = stationary_distribution(vsm1, P_matice) print(f'Stacionární rozdělení (𝜋), a zkouška (𝜋*P)') print(np.concatenate((pi.reshape(-1, 1), (pi.reshape(1, -1) @ P_matice).T), axis=1))

Stacionární rozdělení (𝜋), a zkouška (𝜋*P)
[[0.18148114 0.18148114]
 [0.06850115 0.06850115]
 [0.0110864  0.0110864 ]
 [0.02404211 0.02404211]
 [0.03419535 0.03419535]
 [0.10936119 0.10936119]
 [0.02056555 0.02056555]
 [0.01509768 0.01509768]
 [0.04309803 0.04309803]
 [0.06289474 0.06289474]
 [0.00133551 0.00133551]
 [0.00267245 0.00267245]
 [0.03005804 0.03005804]
 [0.02604366 0.02604366]
 [0.05784417 0.05784417]
 [0.05607889 0.05607889]
 [0.01442454 0.01442454]
 [0.00120198 0.00120198]
 [0.04955565 0.04955565]
 [0.05168766 0.05168766]
 [0.06878944 0.06878944]
 [0.02002891 0.02002891]
 [0.0122929  0.0122929 ]
 [0.01468833 0.01468833]
 [0.00146884 0.00146884]
 [0.02097156 0.02097156]
 [0.00053412 0.00053412]]

3. Porovnejte rozdělení znaků druhého textu se stacionárním rozdělením 𝜋 , tj. na hladině významnosti 5 % otestujte hypotézu, že rozdělení znaků druhého textu se rovná rozdělení 𝜋 z předchozího bodu.

n = len(vsm2.text) konv_tabulka = [] for i in range(len(pi)): konv_tabulka.append([vsm1.letters[i]]) try: konv_tabulka[-1].append(np.round(n*pi[i], 3)) except IndexError: konv_tabulka[-1].append(0) try: konv_tabulka[-1].append(vsm2.letters_abs_freq[i]) except IndexError: konv_tabulka[-1].append(0) konv_tabulka

pi_casted = pi*n t2_letters_counts_casted = vsm2.letters_abs_freq.copy() while np.min(pi_casted) < 5: order = np.argsort(pi_casted) pi_casted[order[1]] += pi_casted[order[0]] t2_letters_counts_casted[order[1]] += t2_letters_counts_casted[order[0]] pi_casted = pi_casted[~np.isin(np.arange(pi_casted.shape[0]), order[0])] t2_letters_counts_casted = t2_letters_counts_casted[~np.isin(np.arange(t2_letters_counts_casted.shape[0]), order[0])] konv_tabulka = np.zeros(( np.max((len(pi_casted), len(t2_letters_counts_casted))), 2)) for i in range(konv_tabulka.shape[0]): try: konv_tabulka[i, 0] = pi_casted[i] except IndexError: pass try: konv_tabulka[i, 1] = t2_letters_counts_casted[i] except IndexError: pass konv_tabulka

plt.bar(np.arange(konv_tabulka.shape[0]), konv_tabulka[:, 0], alpha=0.5, label = 'Stacionární rozdělení pro 1. text') plt.bar(np.arange(konv_tabulka.shape[0]), konv_tabulka[:, 1], alpha=0.5, label = 'Absolutní četnosti znaků ve 2. textu') plt.title("Porovnání výběru s teorií") plt.ylabel('Hodnoty') plt.xlabel("Číslo znaků") plt.legend()

Spočítáme testovou statistiku.

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

Testová statistika: 127.91679723545695

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:" ,pi_casted.shape[0] - 1)

Počet stupňů volnosti: 25

alpha = 0.05 krit_hodnota= stats.chi2.ppf(1-alpha, pi_casted.shape[0] - 1) print("Kritická hodnota:" ,krit_hodnota)

Kritická hodnota: 37.65248413348277

Pokud Testová statistika = 127.917 a Kritická hodnota= 37.652

127.917> 37.652=>

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)

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.chi2.html

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

p-hodnota = 8.692973953159999e-16

.css-ytumd6{-webkit-text-decoration:none;text-decoration:none;}Data a parametry

Nastavení

Výsledky pro první text

Výsledky pro druhý text

1. Za předpokladu výše odhadněte matici přechodu markovského řetězce pro první text. Pro odhad matice přechodu vizte přednášku 17. Odhadnuté pravděpodobnosti přechodu vhodně graficky znázorněte, např. použitím heatmapy.

2. Na základě matice z předchozího bodu najděte stacionární rozdělení 𝜋 tohoto řetězce pro první text.

Výsledky pro první text

3. Porovnejte rozdělení znaků druhého textu se stacionárním rozdělením 𝜋 , tj. na hladině významnosti 5 % otestujte hypotézu, že rozdělení znaků druhého textu se rovná rozdělení 𝜋 z předchozího bodu.

Určení p-hodnoty testu

Data a parametry