I am in the process of learning how to do some basic programming in Mathematica, and so, with my very limited abilities, I attempted to write a simulation in Mathematica today. And I thought the birthday problem seems like a good thing to model.

The problem itself is often phrased something along the lines of:

In a party of N people, how many people we should select, so that there is 50% chance that at least two people shares the same birthday?

We are assuming that there isn’t leap year (no one is born on 29 Feb) and each birthday is equally probable etc.

Few weeks back, I wrote a MATLAB program, and it was basic enough so that when I input a particular *k* value, it would give me a probability out of 100,000 experiments.

============

count=0; times=1e5; k=23; for i=1:times x = round(365*rand(k,1)); if length(x) ~= length(unique(x)) count=count+1; end count; end prob=count/times

============

Just as probability theory predicted, we expect selecting 23 people would be enough to give 50%.

http://en.wikipedia.org/wiki/Birthday_problem

So, that was MATLAB, but I desire a plot of points, which I could vary the value of *k, *something like this:

http://en.wikipedia.org/wiki/File:Birthday_Paradox.svg

But I ran into some problems, because I am not particular familiar with how to manipulate vectors and functions in Mathematica.

==========

f[k_] := ( count = 0; repeat = 10000; Do[ a = RandomInteger[{1, 365}, k]; b = Union[a]; If[Length[a] != Length[b], count++ ] , {i, 1, repeat} ]; prob = N[count/repeat]; Return[prob]; ) list = Table[f[i], {i, 1, 365}]; ListPlot[list] TimeUsed[]

- The whole thing took about 75 seconds to complete. Which is much longer than I expected to be. Starting from 100, the values are so close to probability of 1, it is taking quite a long time to run the simulation. With some playing around, I found most of the time were spent on the “Union” function.
- My original intent was to create a function that list out each probability for each k value. I planned to do something like (as if I am using MATLAB)
list=(list, prob);

HoweverI never figure out how to do this concatenation. (Probably it is for the better, this repeated concatenation might further increase computation time)

- I managed to skip problem 2 with the self defined function f[k_]. However, it took me quite a while to figure out how to use Return[], which gave me exactly what I desire from the function f[k_].

So, the next thing I should do is to consider how to improve the efficiency of the program by modify the “Union” function into something faster.