% Metody teoretyczne

clear
N=1; T=2; % falsz, prawda


% ----------------- zadawanie wartosci prawdopodobienstw

% wlamanie
ps(T)=0.3;
ps(N)=1-ps(T);

% alarm przy warunku wlamania i uderzenia pioruna
pwrs(T,N)=0.6;
pwrs(T,T)=0.0;
pwts(T,N)=0.5;
pwts(T,T)=0.8;
pwrs(N,N)=1-pwrs(T,N);
pwrs(N,T)=1-pwrs(T,T);
pwts(N,T)=1-pwts(T,T);
pwts(N,N)=1-pwts(T,N);


pwhrt(T,N,N)=0.3;
pwhrt(T,N,T)=0.5;
pwhrt(T,T,N)=0.4;
pwhrt(T,T,T)=0.9;
pwhrt(N,N,N)=1-pwhrt(T,N,N);
pwhrt(N,N,T)=1-pwhrt(T,N,T);
pwhrt(N,T,N)=1-pwhrt(T,T,N);
pwhrt(N,T,T)=1-pwhrt(T,T,T);


% ----------------- tablica prawdopodobienstw lacznych

for S=1:2
    for R=1:2
        for Ti=1:2
            for H=1:2
                p(S,R,Ti,H)=pwhrt(H,R,Ti)*pwrs(R,S)*pwts(Ti,S)*ps(S)
            end
        end
    end
end

% Sprawdzenie czy suma wszystkich prawdopodobienstw jest rowna 1
sum(sum(sum(sum(sum(p)))))


% ------------- obliczenie wybranych p-w z rozkładu łącznego
%p(t\h,~s)
% Obliczenie P(A|B,S,U,W)=P(A,B,S,U,W)/P(B,S,U,W)

% Obliczenie P(A|S,B)=P(A,S,B)/P(S,B)
%pwasb=sum(sum(sum(p(2,2,2,:,:))))/sum(sum(sum(sum(p(:,2,2,:,:)))))
%pwnasb=sum(sum(sum(p(1,2,2,:,:))))/sum(sum(sum(sum(p(:,2,2,:,:)))))


pths=sum(sum(p(1,:,2,2)))/sum(sum(sum(p(1,:,:,2))))

pths=sum(sum(sum(p(1,:,2,2))))/sum(sum(p(1,:,2,:)))
% -------------------------------------------------------

% Metody Monte Carlo szacowania wybranych prawdopodobienstw

K=10000; % jednorazowa porcja taktow
srp_as=0; % poczatkowy wynik P(A|S)
srp_a=0; % poczatkowy wynik P(A)

close
hold on
plot([1 1000],[pwas pwas], 'b:')
plot([1 1000],[pa pa],'r:')
axis([0 1000 0 0.1]) % zadanie skali rysunku

% pętla w której będziemy poprawiać naszą estymatę, iteracyjnie licząc średnie p-wo
for j=1:1000

    % losowanie realizacji zgodnie z określonymi w sieci prawdopodobieństwami
    w=rand(1,K)<pw(T);
    u=rand(1,K)<pu(T);

    aNN=rand(1,K)<pwawu(T,N,N);
    aNT=rand(1,K)<pwawu(T,N,T);
    aTN=rand(1,K)<pwawu(T,T,N);
    aTT=rand(1,K)<pwawu(T,T,T);
    a=(w&u&aTT)|(w&~u&aTN)|(~w&u&aNT)|(~w&~u&aNN);

    sN=rand(1,K)<pwsa(T,N);
    sT=rand(1,K)<pwsa(T,T);
    s=(a&sT)|(~a&sN);

    bN=rand(1,K)<pwba(T,N);
    bT=rand(1,K)<pwba(T,T);
    b=(a&bT)|(~a&bN);


    % Obliczenie P(A|S)
    pwas_mc = sum(a&s) / sum(s); % obliczenie czestosci dla biezacej porcji taktow

    % obliczenie czestosci dla wszystkich dotychczasowych taktow
    % (znany wzor na iteracyjne liczenie średniej)
    srp_as = srp_as + (pwas_mc - srp_as) / j;



    % Obliczenie P(A)
    pa_mc = sum(a)/K;
    srp_a = srp_a + (pa_mc - srp_a) / j;




    % aktualizacja wykresu
    historia_sr_as(j) = srp_as;
    plot([1:j], historia_sr_as, 'b');

    historia_sr_a(j) = srp_a;
    plot([1:j], historia_sr_a, 'r')

    pause(0.001)

end