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)
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]}")
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]}")
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))
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)
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)
alpha = 0.05
krit_hodnota= stats.chi2.ppf(1-alpha, pi_casted.shape[0] - 1)
print("Kritická hodnota:" ,krit_hodnota)
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)
print(f"p-hodnota = {scipy.stats.chi2.sf(test_stat, pi_casted.shape[0] - 1)}")