clear all; close all; clc

%% načtení fasta souboru
[header seq] = fastaread('X64775.fasta'); % oryza sativa s mikrosatelity 
% [header seq] = fastaread('AL133493_human_chr21.fasta'); % human chr21 s minisatelity
% [header seq] = fastaread('randomseq_repetice_z_oryzy.fasta'); % náhodná sekvence s mikrosatelity z O. sativa
% [header seq] = fastaread('neisseria gonorrhoeae clone pNG273.fasta'); % rozptýlené repetice bakterie N. gonorrhoeae 
% [header seq] = fastaread('Rabbit short interspersed C repeat.fasta'); % rozptýlené repetice typu SINE (13bp)
% [header seq] = fastaread('Caenorhabditis elegans Cosmid F56F11.fasta'); % gen F56F11 C. elegans s pěti exony
% [header seq] = fastaread('NT_011519.9 human chr22.fasta'); % CpG ostrůvky v chr22
% [header seq] = fastaread('Homo sapiens chromosome 21.fasta'); % CpG ostrůvky v chr21
% [header seq] = fastaread('Mycoplasma pneumoniae.fasta'); % část genomu 1-100000bp
% [header seq] = fastaread('C. carsonella ruddii.fasta'); % celý genom 160000bp
% [header seq] = fastaread('umele_rozptylene_repetice.fasta'); % umělé rozptýlené repetice

% seq = randseq(60000); % náhodná sekvence o délce 60 kbp

%% numerická reprezentace
L = length(seq); % délka sekvence

num_repr = zeros(1,L); 
num_repr(seq =='A') = 0.1260; % metoda EIIP (Electron-Ion Interaction Potential)
num_repr(seq =='T') = 0.1335;
num_repr(seq =='G') = 0.0806;
num_repr(seq =='C') = 0.1340;

%% výpočet kritérií pro výběr hodnoty řádu modelu
q = 70; % max. hodnota řádu modelu - nastavení odhadem

a = zeros(1,q); % vektor parametrů AR modelu
FPE = zeros(1,q); 
AIC = zeros(1,q); 
BIC = zeros(1,q); 
HQC = zeros(1,q); 
    
for i = 1:q
    [a,Vp] = arburg(num_repr,i); % odhad parametrů 'a' AR modelu Burgovou metodou
    % Vp...hodnota ztrátové fce 'Vp' = normalizovaný součet druhých mocnin náhodných chyb e(n)- rozdíl mezi odhadovaným a vstupním signálem
        
    FPE(i) = Vp*((L+i)/(L-i)); % finální predikce chyby
    AIC(i) = log(Vp)+((2*i)/L); % Akaikeho informační kritérium
    BIC(i) = log(Vp)+((i*log(L))/L); % Bayesovo informační kritérium
    HQC(i) = log(Vp)+((2*i*log(log(L)))/L); % Hannan-Quinnovo kritérium
    
    % normalizace osy y na hodnoty 0-1:
    MaxFPE = max(max(FPE)); MinFPE = min(min(FPE));
    nFPE = FPE - MinFPE; nnFPE = nFPE/(MaxFPE-MinFPE);
    
    MaxAIC = max(max(AIC)); MinAIC = min(min(AIC));
    nAIC = AIC - MinAIC; nnAIC = nAIC/(MaxAIC-MinAIC);
    
    MaxBIC = max(max(BIC)); MinBIC = min(min(BIC));
    nBIC = BIC - MinBIC; nnBIC = nBIC/(MaxBIC-MinBIC);
    
    MaxHQC = max(max(HQC)); MinHQC = min(min(HQC));
    nHQC = HQC - MinHQC; nnHQC = nHQC/(MaxHQC-MinHQC);
end

%% vykreslení kriteriálních fcí
plot(nnFPE,'-r+'); hold on
plot(nnAIC,'-go'); hold on
plot(nnBIC,'-b.'); hold on
plot(nnHQC,'-c*'); hold off
title('kritéria pro výběr hodnoty řádu modelu')
xlabel('řád')

legend('FPE','AIC','BIC','HQC','Location','Best')

