from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import scipy.stats as st
import scipy.linalg as la
from tabulate import tabulate
import string
import copy
# Funkce pro nacteni souboru
def load_data(path):
data=""
with open(path) as f:
next(f) # Preskoceni prvniho radku
for line in f:
data += line
return data
# Funkce pro vypocet odhadu pravdepodobnosti a cetnosti znaku v souboru
def symbol_probability(data):
symbols={}
# Pocet kazdeho znaku abecedy (a mezery) nejprve nastavime na 0
# string.ascii_lowercase = 'abcdefghijklmnopqrstuvwxyz'
for s in string.ascii_lowercase:
symbols[s] = 0
symbols[' '] = 0
# Zjisteni poctu pro kazdy znak v datech
for s in data:
if s in symbols:
symbols[s] += 1
# Prevedeni do pravdepodobnosti
control_sum=0
p_symbols={}
for s in symbols:
p_symbols[s] = symbols[s]/len(data)
control_sum += p_symbols[s]
return p_symbols, symbols, control_sum
# Funkce pro graf
def plot(data, title, x, y):
plt.bar(range(len(data)), list(data.values()), align='center')
plt.xticks(range(len(data)), list(data.keys()))
plt.title(title)
plt.xlabel(x)
plt.ylabel(y)
plt.show()
Z obou datových souborů načtěte texty k analýze. Pro každý text zvlášť zjistěte absolutní četnosti jednotlivých znaků (symbolů včetně mezery), které se v textech vyskytují. Dále předpokládejme, že první text je vygenerován z homogenního markovského řetězce s diskrétním časem.
# Nacteni souboru pomoci predchozi funkce
path="source/"
file1="011.txt"
file2="018.txt"
data1=load_data(path+file1)
data2=load_data(path+file2)
print(data1)
print(data2)
# Pravdepodobnosti, cetnosti a kontrolni soucet pro symboly v souboru č.1 a souboru č.2
p1, symbols1, s1=symbol_probability(data1)
p2, symbols2, s2=symbol_probability(data2)
headers=["symbol", "četnost"]
table1=pd.DataFrame(symbols1.items(), columns=headers).to_string(index=False)
table2=pd.DataFrame(symbols2.items(), columns=headers).to_string(index=False)
print("Prvni text:")
print(table1)
plot(symbols1, "Absolutní četnosti znaků v prvním textu", "symbol", "četnost v textu")
print("Druhy text:")
print(table2)
plot(symbols2, "Absolutní četnosti znaků v druhým textu", "symbol", "četnost v textu")
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.
def P_matrix(data):
s_int = {}
n = len(string.ascii_lowercase) + 1
P = np.zeros([n, n])
number = 0
for s in string.ascii_lowercase:
s_int[s] = number
number += 1
s_int[' '] = number
for i in range(0, len(data)):
if i == len(data) - 1:
break
P[s_int[data[i]], s_int[data[i+1]]] += 1
return P/P.sum(axis=1, keepdims=True)
symbols_list = []
for s in string.ascii_lowercase:
symbols_list.append(s)
symbols_list.append(' ')
P = P_matrix(data1)
fig, ax = plt.subplots(dpi=170)
im = ax.imshow(P, cmap='gray', interpolation='nearest')
ax.set_xticks(np.arange(len(symbols_list)), labels=symbols_list)
ax.set_yticks(np.arange(len(symbols_list)), labels=symbols_list)
plt.show()
2. Na základě matice z předchozího bodu najděte stacionární rozdělení π tohoto řetězce pro první text.
def StacDist(P):
W = np.transpose(P-np.eye(P.shape[0]))
pi = la.null_space(W)
pi = np.transpose(pi/sum(pi))
return pi
pi = StacDist(P)
headers=["pi", "hodnota"]
table3=pd.DataFrame(dict(zip(list(range(0, 27)), pi[0])).items(), columns=headers).to_string(index=False)
print(table3)
# Soucet prvku pi je 1
print(sum(pi[0]))
plot(dict(zip(symbols_list, pi[0])), "[Pi] Stacionární rozdělení", "symbol", "hodnota")
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.
table_list = []
for i in range(0, len(pi[0])):
table_list.append([symbols_list[i], len(data1)*pi[0][i], len(data2)*p2[symbols_list[i]]])
headers=["Nij","Pi", "Sigma"]
table=pd.DataFrame(table_list, columns=headers).to_string(index=False)
print(table)
for i in range(0, 3):
table_list[16][i] += table_list[23][i]
table_list[16][i] += table_list[-2][i]
table_list.pop(23)
table_list.pop(-2)
new_table=pd.DataFrame(table_list, columns=headers).to_string(index=False)
print(new_table)
Ni = []
Nj = []
for el in table_list:
Ni.append(el[1])
Nj.append(el[2])
Nij = np.matrix([Nj,Ni])
# Odhady marginal
n = np.sum(Nij)
pi_ = np.sum(Nij, axis = 1)/n
p_j = np.sum(Nij, axis = 0)/n
print("Marginaly")
print("pi_ =\n", pi_)
print("p_j =\n", p_j)
# Teoreticke cestnosti
pipj = np.matmul(pi_,p_j)
npipj = n * pipj
print("Teoreticke cetnosti npipj")
print(npipj)
alpha = 0.05
# Stupen volnosti pomoci rozmeru kontingencni tabulky
df = (np.size(Nij,axis =0) - 1)*(np.size(Nij,axis =1) - 1)
print(f"stupeň volnosti: {df}")
# Kriticka hodnota s alpha a df
chi2 = st.chi2.isf(alpha,df)
print(f"kriticka hodnota: {chi2}")
# Testova statistika
Chi2 = np.sum(np.square(Nij - npipj)/npipj)
p = st.chi2.sf(Chi2,df)
print(f"testová statistika: {Chi2}")
print(f"p-hodnota: {p}")
if Chi2 > chi2:
print("Zamitame H0")
else:
print("Nezamitame H0")