function [ probabilities ] = NCoincidentBirthdaysProbability( maxBirthdays, numDays, numTrials, maxGroupSize )
probabilities = zeros(maxBirthdays-1, maxGroupSize);
for nBirthdays = 2:maxBirthdays
    for nPeople = 1:maxGroupSize
        prob = zeros(1, numTrials);
        for i = 1:numTrials
            birthdays = ceil(numDays*rand(1, nPeople));
            anyRepeated = length(unique(birthdays) < length(birthdays));
            if anyRepeated
                uniqueValues = unique(birthdays);
                for j = 1:numel(uniqueValues)
                    IDX = find(birthdays == uniqueValues(j));
                    if numel(IDX) >= nBirthdays
                        prob(i) = 1;
                    end
                end
            end
        end
        probabilities(nBirthdays-1, nPeople) = sum(prob)/numTrials;
    end
end
end