%
%     Bernoulli.M
%
% Analyse des Bernoulliprozesses
%
% Jan Lunze
% für 3. Auflage: 5.1.2017
% für 4. Auflage: 1.8.2025
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% 
clear
close all

%%  Geometrische Verteilung
%  Prob(X=n) = (1-p) p^n
figure(1)
q=0.3;
N1=5;
NAchse=0:N1;
P=[];
k1=0;
for k2=NAchse
    k1=k1+1;
    P(k1)=q*(1-q)^k2;
end
subplot(2,1,1)
stairs(NAchse, P, 'b', 'LineWidth',1);
hold on
plot(NAchse, P, 'bo');
latexxlabel('$$m$$');
xticks([0 1 2 3 4 5 6]);
latexylabel('$${\rm Prob}(X=m)$$');
rotateY;
yticks([0 0.3]);
axis([0 N1 0 1.1*q]);
hold off
    

%%  Binomialverteilung
%   B(z|p,k) = (k/z) p^z (1-p)^{k-z}
%
figure(2)
p=0.7;
N1=10;
NAchse=0:N1;
% z=0
z=0;
P=[];
k1=0;
for k2=NAchse      % Zeitpunkte k2=0,...,N1
    k1=k1+1;
    P(k1)=factorial(k2)/factorial(z)/factorial(k2-z)*p^z*(1-p)^(k2-z);
end
subplot(3,1,1)
stairs(NAchse, P, 'b', 'LineWidth',1);
hold on
plot(NAchse, P, 'bo');
%latexxlabel('$$k$$');
%xticks([0 1 2 3 4 5 6]);
xleer;
latexylabel('$$B(0\,|\,0,\!7, k)$$');
rotateY;
yticks([0 1]);
xticksempty([0 5 10]);
axis([0 N1 -0.2 1.1]);
hold off
% z=2
z=2;
P=[];
for k1=1:z        % Zeitpunkte k=0,...,z-1; dargestellt durch k1=1,...,z
    P(k1)=0;
end
for k2=z:N1    % Zeitpunkte k2=z,...,N1; 
    k1=k1+1;
    P(k1)=factorial(k2)/factorial(z)/factorial(k2-z)*p^z*(1-p)^(k2-z);
end
%
subplot(3,1,2)
stairs(NAchse, P, 'b', 'LineWidth',1);
hold on
plot(NAchse, P, 'bo');
%latexxlabel('$$k$$');
%xticks([0 1 2 3 4 5 6]);
xleer;
latexylabel('$$B(2\,|\,0,\!7, k)$$');
rotateY;
yticks([0  1]);
xticksempty([0 5 10]);
axis([0 N1 -0.2 1.1]);
hold off
% z=3
z=3;
P=[];
for k1=1:z        % Zeitpunkte k=0,...,z-1; dargestellt durch k1=1,...,z
    P(k1)=0;
end
for k2=z:N1    % Zeitpunkte k2=z,...,N1; 
    k1=k1+1;
    P(k1)=factorial(k2)/factorial(z)/factorial(k2-z)*p^z*(1-p)^(k2-z);
end
%
subplot(3,1,3)
stairs(NAchse, P, 'b', 'LineWidth',1);
hold on
plot(NAchse, P, 'bo');
latexxlabel('$$k$$');
xticks([0 5 10]);
latexylabel('$$B(3\,|\,0,\!7, k)$$');
rotateY;
yticks([0 1]);
axis([0 N1 -0.2 1.1]);
hold off

%%   Vergleich von Erwartungswert und einer Realisierung
%
figure(3)
subplot(3,1,[1 2])
%  Berechnung des erwarteten Verhaltens
p=0.7;
N1=20;
NAchse=0:N1;
plot([0 N1],[0 N1*p],'b--','LineWidth',1.5);  % Gerade für den Erwartungswert
hold on
%  Berechnung einer Realisierung
load('RandBernoulli');       % Anfangswert des Zufallsgenerators
rng(RandBernoulli);
R=rand([1 N1+1]);   %Zufallszahlen im Intervall (0,1)
Z=[];
Z(1)=0;
for k1=1:N1
    if R(k1) < p        % wenn Zufallszahl größer als p: Zustandswechsel
        Z(k1+1)=Z(k1)+1;
    else
        Z(k1+1)=Z(k1);
    end
end
stairs(NAchse, Z, 'b', 'LineWidth',1.5);
plot(NAchse, Z, 'bo');
axis([0 N1 0 15]);
latextitle('Bernoulliprozess ($$p=0.7$$)');
latexxlabel('$$k$$');
xticks([0 5 10 15 20 25]);
latexylabel('$$z$$');
rotateY;
yticks([0 5 10 15]);
hold off

%%
epsfigc('Bernoulli');
