matlaby siem zbiegają

% Metody teoretyczne

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

%zadawanie wartosci prawdopodobienstw

%Pr - błąd programisty
pPr(T) = 0.70;
pPr(N) = 1-pPr(T);

%N - skok napięcia
pN(T) = 0.07;
pN(N) = 1-pN(T);

%B - błąd w systemie
pBpr(T, T) = 0.13; %blad programisty -> blad w systemie
pBpr(N, T) = 0.05;
pBpr(N, N) = 1 - pBpr(N, T);
pBpr(T, N) = 1 - pBpr(T,T);

% bład na dysku przy skoku napiecia i bledzie programisty
pDnb(T,T,T) = 0.80;
pDnb(T,T,N)= 1-pDnb(T,T,T);
pDnb(N,T,T)= 0.15;
pDnb(N,T,N) = 1-pDnb(N,T,T);
pDnb(T,N,T) = 0.45;
pDnb(T,N,N) = 1-pDnb(T,N,T);
pDnb(N,N,T) = 0.06;
pDnb(N,N,N) = 1-pDnb(N,N,T);

% utrata danych przy bledzie na dysku i bledzie w systemie
pUdb(T,T,T) = 0.5;
pUdb(T,T,N) = 1 - pUdb(T,T,T);
pUdb(N,T,T) = 0.3;
pUdb(N,T,N) = 1 - pUdb(N,T,T);
pUdb(T,N,T) = 0.2;
pUdb(T,N,N) = 1 - pUdb(T,N,T);
pUdb(N,N,T) = 0.05;
pUdb(N,N,N) = 1 - pUdb(N,N,T);

% tablica prawdopodobienstw lacznych

for Pr=1:2
  for B=1:2
    for Ni=1:2
      for D=1:2
        for U=1:2
          p(Pr,B,Ni,D,U)=pPr(Pr)*pN(Ni)*pBpr(Pr,B)*pDnb(Ni,B,D)*pUdb(D,B,U);
        end
      end
    end
  end
end

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


% Na podstawie tablicy prawdopodopienstw lacznych mozna obliczyc dowolne prawdopodobienstwo, np.P(A|S)=P(A,S)/P(S)

%pwas=sum(sum(sum(p(2,:,2,:,:))))/sum(sum(sum(sum(p(:,:,2,:,:)))))
%pwnas=sum(sum(sum(p(1,:,2,:,:))))/sum(sum(sum(sum(p(:,:,2,:,:)))))
pD=sum(sum(sum(sum(p(:,:,:,2,:)))));
pPrd=sum(sum(sum(p(2,:,:,2,:))))/sum(sum(sum(sum(p(2,:,:,:,:)))));
pD
pPrd
% -------------------------------------------------------
% Metody Monte Carlo

K=10000; % jednorazowa porcja taktow
close
hold on                             % polecenie zostawia wykres, nastepne beda rysowane na tym samym
plot([1 10],[pPrd pPrd])          % wykres P(A|S) wyznaczonego teoretycznie
axis([0 10 0 1])                % zadanie skali rysunku

sredniaD(1)=0;
sredniaPrD(1)=0;
for j=1:1000
    %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);

    pr=rand(1,K)<pPr(T);
    n=rand(1,K)<pN(T);
    bT=rand(1,K)<pBpr(T,T);
    bN=rand(1,K)<pBpr(N,T);
    bb=(pr&bT)|(~pr&bN);

    dNN=rand(1,K)<pDnb(N,N,T);
    dNT=rand(1,K)<pDnb(N,T,T);
    dTN=rand(1,K)<pDnb(T,N,T);
    dTT=rand(1,K)<pDnb(T,T,T);
    dd=(n&bb&dTT)|(~n&~bb&dNN)|(n&~bb&dTN)|(~n&bb&dNT);
    % ......itd.............

    % obliczenie czestosci w j-tym kroku na podstawie odpowiednich tablic jedynek, np. P(A|S)= sum(a&s)/sum(s)
    PD=sum(dd)/length(dd);
    PPrD=(sum(dd&pr)/length(dd))/(sum(pr)/length(pr));

    sredniaPrD(j+1) = sredniaPrD(j)+ (PPrD-sredniaPrD(j))/j;
    sredniaD(j+1) = sredniaD(j) + (PD-sredniaD(j))/j;
    % usrednienie czestosci dla wszystkich dotychczasowych krokow
plot([1:j+1], sredniaPrD);
    % narysowanie wykresu: plot([1:j], tablica_srednich_czestosci)
    %plot([1:j+1], sredniaD);
    %pause

end
sredniaPrD(1001)
sredniaD(1001)