%% Geburtstags Problem: 
% Wie hoch ist die Wahrscheinlichkeit, dass bei n Personen in einem Raum
% mindestens zwei davon am gleichen Tag Geburtstag haben?

% Anzahl Tage im Jahr
d = 365;

% Was soll die groesste Gruppe sein, welche betrachtet wird
max_people = 70;

% Die Berechnung wird fuer alle Gruppen-Groessen von 1 .. max_people
% gemacht
N = 1:max_people;

%% Berechnung durch Betrachtung der Wahrscheinlichkeiten dass niemand den gleichen Geburtstag hat

hold off

% Initialisierung aller Wahrscheinlichkeiten auf 1 (neutrales Element der
% Multiplikation)
P = ones(1,max_people);

% Klein n geht also von 1 .. max_people.
% Fuer jedes n berechnen wir die gesuchte Wahrscheinlichkeit.
for n = N
    % Berechnung der Reihe 365/365 * 364/365 * ... * (365 - n + 1)/365
    for i = 1 : n
        P(n) = P(n) * (d-i+1) / d;
    end
    % Die Wahrscheinlichkeit dass mindestens zwei gleich sind ist 1 - die berechnete Wahrscheinlichkeit
    P(n) = 1 - P(n);
end

plot(P)
grid on
xlabel('Anzahl Personen','FontSize',20)
ylabel('Wahrscheinlichkeit','FontSize',20)

hold on

%% Simulation zur Ueberpruefung der Formel

at_least = 2;
experiments = 1000;

P_sim = zeros(1,max_people);

% Simuliere Gruppen von 1 .. max_people.
for n = N
    
    % Zaehlt die Experimente in denen mindestens 'at_least' personen am
    % gleichen Tag Geburtstag haben
    successful = 0;
    
    for i = 1:experiments
        
        % Kumuliere die Geburtstage pro Tag im Jahr, initial 0.
        BirthdayCount = zeros(1, d);
        
        % Waehle 'n' zufaellige Geburtstage
        Birthdays = randi([1, d],1,n);
        
        % Fuer jeden Geburtstag, erhoehe den Zaehler des ensprechenden
        % Tages um 1
        for birthday = Birthdays
            BirthdayCount(birthday) = BirthdayCount(birthday) + 1;
        end
        
        % Inkrementiere 'successful' falls in diesem experiment mindestens
        % 'at_least' personen am selben Tag Geburtstag haben.
        maximum = max(BirthdayCount);
        if maximum >= at_least
            successful = successful + 1;
        end
    end
    
    % Die gemessene Wahrscheinlichkeit
    P_sim(n) = successful/experiments;
end

% Glaettung der Punkte. Wissenschaftlich wahrscheinlich  nicht gut.
%P_sim = smooth(P_sim,7);

plot(P_sim, ':.','MarkerSize',10)
axis([0 max_people 0 1])
grid on
xlabel('Anzahl Personen','FontSize',20)
ylabel('Wahrscheinlichkeit','FontSize',20)

% Zeige an, wie viele Random number gebraucht wurden fuer diese Simulation.
disp('Number of generated random numbers: ')
disp(max_people * (max_people-1)/2*experiments);

%% Clear graphics

hold off