# Probability: proportion of 1000 random lists for x that contain the same nrs

I want to estimate, for groups of x persons, the probability that at least two of them have the same birthday. I want to repeat this process 1000 times. So I'm looking for the number of 1000 random birthday lists that contain the same birthday. (in other words, I need the probability)

Let' say a function which calculates a random birthdays is ($x$ can be from $x=1$ to $x=365$)

birthdays[x_] := RandomInteger[{1, 365}, x]


is my answer. Take x=25 and determine whether birthdays[x] contain equal numbers.

bday1 = birthdays[25]

Length[bday1]

ubday1 = Union[bday1]


Length[ubday1] is my answer.

This is my answer so far; and I don't know whether its correct or not.

birthprob[x_] := Count[Table[RandomInteger[{1, 365}, x], {1000}],


Could someone help me?

-
What do the following shorthands mean: "bdays" and "nrs"? I think I know what they are, but it is best if you post complete words and not shorthand notation. – rcollyer Oct 10 '12 at 15:32
sorry, it means birthdays and numbers :) as you probably expected – Jaimy Oct 10 '12 at 15:33
You question is not clear, put some example. – Murta Oct 10 '12 at 15:39
Yes, I did expect that is what it meant, but while shorthand works for instant messaging, it is not, in general, good for overall communication. So, if you would update your question using this link without the shorthand, with as complete sentences as you're capable of (if your a non-native speaker, we can edit it into shape), and with examples of what you tried. – rcollyer Oct 10 '12 at 15:42
You might be interested in the functions Count[], Tally[], and Gather[], as well as Histogram[]. – J. M. Oct 10 '12 at 16:09

This is not particularly different to any of the excellent answers which I have all upvoted. It did, however, prompt a limited comparison (benefiting from the workaround for bug in BenchmarkPlot here. Testing against @Mr.Wizard birthprob2 (just because I liked it):

fun[n_, m_] := With[{rv = RandomInteger[{1, 365}, {m, n}]},
Total@(Boole[Length@Union[#] < Length[#]] & /@ rv)/m]
birthprob2[x_Integer?Positive, n_Integer?Positive] :=
1 - Sum[Abs@Signature@RandomInteger[{1, 365}, x], {n}]/n
Needs["GeneralUtilities"]
Get["GeneralUtilities"];
indexBy[f_][expr_] := expr~indexBy~f
indexBy[asc_Association, f_] := Values[asc]~indexBy~f
p : indexBy[expr_, f_] /; ! AtomQ[expr] ||
Message[indexBy::normal, 1, HoldForm@p] :=
AssociationThread[f /@ expr, expr]
p : indexBy[_, _, __] /;
Message[indexBy::argt, indexBy, Length@Unevaluated@p, 1, 2] := Null
With[{dv :=
DownValues[
GeneralUtilitiesBenchmarkingPackagePrivateaddlabels]},
dv = dv /. SystemIndexBy -> indexBy;]
bp = BenchmarkPlot[{fun[10, #] &, birthprob2[10, #] &, fun[25, #] &,
birthprob2[25, #] &, fun[50, #] &, birthprob2[50, #] &,
fun[200, #] &, birthprob2[200, #] &}, # &, PowerRange[10, 1000000],
"IncludeFits" -> True]


A slight advantage (but all linear) for Union approach.

For completeness:

DiscretePlot[fun[j, 1000], {j, 1, 365}]


-
@Mr.Wizard I just liked your use of Signature and tested against this. This was good fortune as it led me to your question (and answer) re: problem with BenchmarkPlot...so thank you for that post :) – ubpdqn May 21 at 10:23

Your own code is nearly complete. You just need a pattern that checks if a list contains duplicates. I shall use Signature, but you can use any of the methods referenced. I shall also use the syntax of RandomInteger that does away with Table.

birthprob[x_] :=
Count[RandomInteger[{1, 365}, {1000, x}], _?(Signature@# === 0 &)] / 1000

birthprob[23]

0.503


If you try to use this code for much more than 1000 samples you may run out of memory because you are producing a single massive array for all trials. If instead we determine the result for each trial separately we can use a much larger sample size. I shall use Sum, which has a syntax similar to Table. Abs is needed because Signature may return either -1 or 1 in addition to 0.

birthprob2[x_Integer?Positive, n_Integer?Positive] :=
1 - Sum[Abs @ Signature @ RandomInteger[{1, 365}, x], {n}]/n

birthprob2[23, 1*^6] // N  (* 1,000,000 trials *)

0.506157

-

Estimating probabilities via proportions in simulation outcomes is conveniently done by averaging indicators. An "indicator" equals $1$ when an event happens and $0$ when it does not: the average of a list of such values is, mathematically, the proportion of $1$'s in the list.

Your work can be expressed as

birthdays[x_Integer] /; x >= 1 := RandomInteger[{1, 365}, x];
hasDuplicate[x_] := x != Length[Union[birthdays[x]]];


Boole is Mathematica's indicator function and--as one might guess--Mean will average elements of a list. Let's use these, along with your idea to deploy Table for iteration over the simulations, to write a function that finds the proportion of times, out of $n$ independent trials, that $x$ birthdays have a duplicate:

prob[x_Integer, n_Integer] /; n >= 1 := Mean[Table[Boole[hasDuplicate[x]], {n}]]


Apply it:

prob[25, 10^5]


$\frac{11381}{20000}$

Your answer will vary randomly from this one--but it will be close (depending on how large you make $n$).

After sneaking up on such a solution by creating its component parts in steps (birthdays, hasDuplicate, and then prob), some people like to combine it into a single function, which might look something like

prob2[x_Integer, n_Integer] /; x >= 1 && n >= 1 :=
Mean[Boole[x != Length[Union[#]]] & /@ RandomInteger[{1, 365}, {n, x}]];


Although this often makes the solution method less clear, it may improve the computation time and can package it in a form more suitable for further analyses:

SeedRandom[17]; prob[25, 10^6] // N // Timing


{6.411, 0.568418}

SeedRandom[17]; prob2[25, 10^6] // N // Timing


{1.529, 0.568418}

In this case, there is a four-fold speedup. As a simple check that no error was made at this last step, both solutions give the same answer when starting with the same random seed.

-
Wow, i like the short notation needed for a -in my opinion- quite hard problem! Many thanks (everybody who helped me!) – Jaimy Oct 10 '12 at 17:28
While you were writing that comment Jaimy, I was completing an edit that anticipates and reacts to it: there often is a trade-off between clarity and efficiency in Mathematica. I believe a good strategy is to strive for clarity, because that increases the chances of producing correct, maintainable, extendable code, but to be ready when necessary to make the (usually mechanical) transformations that capitalize on some of Mathematica's built-in efficiencies, at the cost of some readability in the code. – whuber Oct 10 '12 at 17:33

I think the full description of the problem you want to solve, is:

Generate 1000 lists, each of which contain x random numbers in the range 1 to 365, and determine how many of them contain duplicate numbers, that is, have the same number appear twice in the list. From the name birthprob I assume that you actually want the (estimated) probability, so you have to divide the obtained number by 1000.

Your solution is incomplete (the second argument of the Count is missing), but it is indeed the initial part of a correct solution. Let me however point out a simplification: Instead of Table[RandomInteger[{1, 365}, x], {1000}] you can simply write RandomInteger[{1, 365}, {1000, x}].

What is missing in your solution is the second element of Count (and, of course, the division by 1000). Count expects as second element a pattern. What you want is a pattern that matches lists with duplicate elements. This pattern is {___, a_, ___, a_, ___}. Therefore the complete solution (with the simplified table generation) reads

birthprob[x_] :=
Count[RandomInteger[{1, 365}, {1000, x}], {___, a_, ___, a_, ___}]/1000


OK, so how does this pattern work? First, ___ matches any sequence of expressions, especially any sequence of numbers, including the empty one. On the other hand, _ matches exactly one expression. Now a_ is a (short-hand) notation to say "I want whatever this pattern matches to be named a" (note that you already used the very same form in the definition of birthprob: by using x_ as the argument, you say "accept any single argument, and name it x). So the full expression reads as "a list which contains an arbitrary number of elements (including none), then an element which shall be named a, then again an arbitrary number of elements (including none), then the single element a again, followed by an arbitrary number of elements". Therefore a list matches this pattern exactly if it contains a duplicate element (which will be matched with a_).

Let's look at this for x up to 100:

ListPlot[Array[birthprob,{100}]]


-
+1 This is an elegant way to translate the statement of the problem directly into code. (I bet it took a while to plot that graph, though. :-) – whuber Oct 10 '12 at 17:39

Assuming that you want the proportion of lists of length x of random integers between 1 and 365 that contain the same element twice or more, you can do something like

birthdays[x_] := RandomInteger[{1, 365}, x]
theAnswer[x_, iters_] := Count[
Count[
Tally[#],
{_, y_} /; y >= 2
] & /@ Table[birthdays[x], {iters}],
z_ /; z > 0
]/iters // N


or

theAnswer2[x_, iters_] :=
Count[Count[Tally[#], {_, y_} /; y >= 2] & /@ Table[birthdays[x], {iters}], z_ /; z > 0]/iters // N


(exactly the same but in one line so it looks less daunting), then eg

Table[{i, theAnswer[i, 1000]}, {i, 1, 100}] //
ListPlot[#, Joined -> True] &


But maybe I misunderstood (or made a mistake in the code).

This works as follows: the

Count[
Tally[#],
{_, y_} /; y >= 2
] &


bit is a function that takes a list of the form returned by Tally[birthdays[x]], that is, a list of lists, each of which has length two, and counts how many have a second element larger than or equal to 2. Then, we map this over a list of iters birthday lists (iters=1000 in your example), obtaining a list of length iters, each element of which corresponds to the number of birthdays occuring more than once in each birthday list. We then count how many elements of that list are larger than 0, ie, how many of the birthday lists had more than 0 coinciding birthdays.

-
Okay i can follow but not reproduce. its my second week with mathematica. and this doesnt give me the answer for the question (this is nice as well), since i have to compute a probability, can you help me with that? – Jaimy Oct 10 '12 at 16:34
what I am plotting is the "experimental estimate for the probability" (vertical axis) versus the length of the list, which you call x (horizontal axis). Basically I'm plotting what is also plotted in the first figure in the wikipedia article. So, theAnswer[25,10000]` (say) gives you the answer (ahem) you are looking for. – acl Oct 10 '12 at 16:44