%Birthday problem simulation method

numPeople = 86;

%Let's ignore leap years
numDays = 356;
numTrials = 100000;
cumTriples = zeros(numTrials, 1);

numPeople = round(numPeople);

for i = 1:numTrials
    birthdays = round(numDays .* rand(numPeople, 1));
    numTriples = sum(hist(birthdays, unique(birthdays)) > 2);
    cumTriples(i) = numTriples;
end

percent = sum(logical(cumTriples)) / numTrials;

% It seems as though 86 is the number of people one would expect to require
% based on this simulation.


% Birthday problem Poisson Approximation

probability = 1 - exp(-(nchoosek(numPeople, 3)) / 365^2);

fprintf('In Simulation 3 in %d people share a birthday with probability %f\n', numPeople, percent);
fprintf('By a Poisson Approximation, 3 in %d people share a birthday with probability %f \n', numPeople, probability);

% It seems as though 84 is the number of people one would expect to require
% if we were going on the Poisson approximation.