echo off; cla; clear; format compact; format short;
'Uloha c. 2'
'Priklad 1: Binomicke rozdeleni pravdepodobnosti'

'a / b = bilych / cernych kulicek'
a = 1
b = 5
'pocet tazeni v ramci pokusu'
s = 16

X = 0:s; % mnozina vsech moznych vysledku pokusu
         % (pocet vytazeni bile kulicky)
p = a/(a+b); % pravdepodobnost vytazeni bile kulicky
q = 1-p;  % pravdepodobnost vytazeni cerne kulicky

% 1) Pouziti statistickeho toolboxu:
    'distribucni funkce rozdeleni pravdepodobnosti'
cdf = binocdf (X, s, p)
    'frekvencni funkce rozdeleni pravdepodobnosti'
pdf = binopdf (X, s, p)

% 2) Vlastni vypocet clenu binomickeho rozvoje:
for i=0:s % pocet vytazeni bile kulicky
  % nize: vzorec pro vypocet pravdepodobnosti vytazeni "i" bilych kulicek pri "s" tazich,
  % pri pravdepodobnosti vytazeni bile kulicky "p"
mypdf (i+1,1) = ( ( factorial(s)/factorial(i)/factorial(s-i) ) * p^(i) * q^(s-i) );
    % cast pro vytvoreni cdf z pdf
  if i == 0
    mycdf (i+1,1) = mypdf (i+1,1);
  else
    mycdf (i+1,1) = mypdf (i+1,1) + mycdf (i,1);
  end
end

% Vypocet smerodatnych odchylek:
pokusu100 = 100;
pokusu1000 = 1000;
for i=1:s+1
  odchylka100 (i,1) = ( pdf(i)*(1-pdf(i))/pokusu100 )^(1/2);
  odchylka1000 (i,1) = ( pdf(i)*(1-pdf(i))/pokusu1000 )^(1/2);
end
for i=1:s+1
  odch_hor100(i) = pdf(i)+odchylka100(i); 
  odch_dol100(i) = pdf(i)-odchylka100(i); 
  odch_hor1000(i) = pdf(i)+odchylka1000(i); 
  odch_dol1000(i) = pdf(i)-odchylka1000(i); 
end

% Vygenerovani 100 pokusu
  % Jednotlive pokusy jsou v radku
  % "1" odpovida vytazeni bile kulicky, "0" cerne
tazeni = (ceil((a+b)*rand(pokusu100,s)));
for i=1:pokusu100
  for j=1:s
    if tazeni(i,j)<=a
      tazeni(i,j)=1;
    else
      tazeni(i,j)=0;
    end
  end
end

% Pravdìpodobnost vytažení "j" bilych kulicek
% Rozhodnutí o ne/pøíznivém pokusu:
for i=1:pokusu100
  for j=1:s+1
    if sum (tazeni(i,:)) == j-1
      priznivych(i,j)=1;
    else
      priznivych(i,j)=0;
    end
  end
end

for j=1:s+1
  pravdepodobnost (j,1) = sum (priznivych(:,j)) / pokusu100;
end

hold on
xlabel('Pocet tahu')
ylabel('Pravdepodobnost')

plot (X, odch_hor1000, ':c')
plot (X, odch_dol1000, ':c')
plot (X, odch_hor100, ':k')
plot (X, odch_dol100, ':k')
plot (X, pdf, '+r')
plot (X, cdf, '-g')
plot (X, pravdepodobnost, 'xb')
hold off

pause
'Priklad 2: Podminena pravdepodobnost'

'a / b = bilich / cernych kulicek'
a = 1
b = 5
'pocet tahu'
s = 5

% Aby nevznikala záporná èísla, nezmeni se pomìr a / b
if a < s | b < s
  a = a*10;
  b = b*10;
end

n = a + b;

% Matice p: radek + 1 ~ tah, sloupec ~ jiz vytazenych bilych
for i=0:s-1 
  for j=0:s-1
    p(i+1,j+1) = ((a-j)/(n-i))* ...
      (factorial(a)/factorial(j)/factorial(abs(a-j)))* ...
      (factorial(b)/factorial(abs(i-j))/factorial(abs(b-(i-j))))/ ...
      (factorial(n)/factorial(i)/factorial(abs(n-i)));
  end
end

'pi - pravdepodobnost v i-tem tahu:'
format rat;
p1 = p(1,1)
p2 = p(2,1) + p(2,2)
p3 = p(3,1) + p(3,2) + p(3,3)
p4 = p(4,1) + p(4,2) + p(4,3) + p(4,4)
p5 = p(5,1) + p(5,2) + p(5,3) + p(5,4) + p(5,5)

format;

% Skript vytvoreny pro vlastni kontrolu (neni soucasti zadani):

% 'Generovani 5 x 1000 tahu:'
% 
% pokusu = 1000;
% 
% vyt = zeros (pokusu,1);
% suma = zeros (s,1);
% opak = zeros (pokusu,s);
% 
% for i=1:pokusu
%   for j=1:s
%     tazeni = [ones(1,a-sum(vyt(i,:))), zeros(1,b-sum(opak(i,:)))];
%     vyt (i,j) = tazeni (1,ceil((n+1-j)*rand));
%     if vyt (i,j) == 1
%       suma (j,1) = suma (j,1) + 1;
%     else
%       opak (i,j) = 1;
%     end
%   end   
% end
% 
% % Vypocet pravdepodobnosti
% gen_p1 = suma(1) / pokusu
% gen_p2 = suma(2) / pokusu
% gen_p3 = suma(3) / pokusu
% gen_p4 = suma(4) / pokusu
% gen_p5 = suma(5) / pokusu
% 
% 'Porovnani pi-gen_pi:'
% format long
% delta1 = p1 - gen_p1
% delta2 = p2 - gen_p2
% delta3 = p3 - gen_p3
% delta4 = p4 - gen_p4
% delta5 = p5 - gen_p5
