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

Question

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?

Answer 1

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.

Answer 2

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}]]

Answer 3

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]
ClearAll[theAnswer];
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.

Answer 4

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

Answer 5

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[
    GeneralUtilities`Benchmarking`PackagePrivate`addlabels]}, 
 dv = dv /. System`IndexBy -> 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}]