%-------------------- cviko 2: HMM  - je to v Matlabu, ale lepsi je si to 
% zkouset a malovat NA PAPIRE S KALKULACKOU ! ---------------------------   
% matlab vyuzit maximalne na vypocet bj[o(t)], to je na papire opravdu dost 
% opruz. 

% pozorovaci sekvence:
T = 5; 
o1=[0.5  1.2 1  -0.5 -1];
o2=[1     3  2  -3   -5];
O=[o1; o2];

% parametry modelu 
N = 4; 
A = [ 0       1       0        0;
      0       0.6     0.4      0;
      0       0       0.7      0.3;
      0       0       0        0    ];
MI2=[1;2]; SIGMA2=[2;2];
MI3=[-1;-2]; SIGMA3=[2; 2]; 
% take tyto hodnoty nacpu do matic, kde prvni a posledni sloupec budou uplne 
% na houby:
MI =    [zeros(2,1) MI2    MI3    zeros(2,1)]; 
SIGMA = [zeros(2,1) SIGMA2 SIGMA3 zeros(2,1)]; 

% jak vypada O
subplot(111);
plot (1:5, O); 

% pocitani vysilacich pravdepodobnosti - radky 1 a 4 budou v te matici 
% uplne na prd (prvni a posledni stav nejsou vysilaci), ale nechme je tam - 
% bude se snaze indexovat !
EMIS=zeros(N,T);
for t=1:T,
  % budu pocitat jen pro vysilaci stavy:
  for j = 2:N-1,
    EMIS(j,t) = normal (O(:,t), MI(:,j), SIGMA(:,j));
  end
end
EMIS

%--- vysledek ---
%EMIS =
%
%    0.0340    0.0349    0.0398    0.0013    0.0001
%    0.0098    0.0010    0.0033    0.0340    0.0129
% zkuste si to spocitat aspon pro 1 vektor a jednu funkci RUCNE !!!!

% stavove sekvence musi mit T+2 prvku, zacinaji 1, konci N, musi respektovat
% matici prechodovych pravdepodobnosti

x1 = [1 2 3 3 3 3 4];
x2 = [1 2 2 3 3 3 4];
x3 = [1 2 2 2 3 3 4];
x4 = [1 2 2 2 2 3 4];

% vic jich neni, dame je vsechny do matice, jejich pocet znacim K
X = [x1; x2; x3; x4]; 
K = 4;

% definuji nejaky vektor pro vsechny P (O,X | M)
POXM = zeros(4,1); 

% a jdeme je spocitat (vzorecek 11)
for k = 1:K
  x = X(k,:); % vybiram stav. sekvenci - pozor, stav sekvence by se mela 
              % indexovat od 0, ale to v matlabu nejde ! 
  P = A(x(1),x(2)); % melo by byt a(x0,x1), ale opakuji, jsme v matlabu
  for t = 1:T,
    xt = x(t+1); xtplus1 = x(t+2); % zase ten matlab
    ot = O(:,t);                   % vybiram vektor o(t)
    bjot = normal (ot, MI(:,xt), SIGMA(:,xt)); %disp(bjot)
    axtxtplus1 = A(xt,xtplus1); %disp(axtxtplus1)
    P = P * bjot * axtxtplus1; 
  end
  % vrazit tam, kam je potreba:
  POXM(k) = P;
end

% vysledek
%POXM =
%   1.0e-08 *
%
%    0.0007
%    0.0211
%    0.2201
%    0.0073

% kontrola treba pro tu druhou stav. sekvenci:
kontrola2 = 1.0 * 0.0340 * 0.6 * 0.0349 * 0.4 * ...
            0.0033 * 0.7 * 0.0340 * 0.7 * 0.0129 * 0.3

%kontrola2 = 6.0592e-11 - zda se ok. 

%%%%%%%%%%% Baum-welch pravdepodobnost: P(O|M) %%%%%%%%%%
POM = sum(POXM)

% 7.1202e-10

%%%%%%%%%%% Viterbiho pravdepodobnost: P*(O|M) %%%%%%%%%%
PstarOM = max(POXM)

% 6.2896e-10 

% ......... a zbytek sami ... muzete vyuzit funkci, ktere jsou na
% www strankach (reestim a viterbi_log)