8
$\begingroup$

[this probably belongs on cross validated ... ]

Working through Mosteller's Fifty Challenging Problems in Probability with Mathematica. I am trying to understand Mosteller's answer for The First Ace problem (#40). The problem goes:

Shuffle an ordinary deck of 52 playing cards containing 4 aces. Then turn up cards from the top until the first ace appears. On the average how many cards are required to produce the first ace?

I solved this problem in Mathematica as as follows:

suite = Flatten@{Range[2, 10], J, Q, K, A}
deck = Flatten@{suite, suite, suite, suite}
dat = Table[First@Flatten@Position[
         RandomSample[deck, Length@deck], A], {n, 1, 100000}];
ListPlot[
    Tally@dat,
    GridLines -> Automatic,
    GridLinesStyle -> Directive[Dotted, Gray],
    Frame -> True
]

enter image description here

Mean@dat // N
10.6152

Which yields Mosteller's answer (by symmetry) of 10.6 cards (48/5 +1).

I was a bit shocked that the mode was 1 card.

Finally, when I look at the CDF the 50% mark is at 8 cards

t = Reverse[SortBy[Tally@dat, #[[2]] &]][[All, 2]]
ListPlot[N[Accumulate@t/100000],
    GridLines -> Automatic,
    GridLinesStyle -> Directive[Dotted, Gray]
]

enter image description here

So my question is, if I played this game for real in a casino, should I go for 8 or 10 cards?

$\endgroup$
3
  • $\begingroup$ I think you're right, this probably belongs on CrossValidated... Cool problem though. $\endgroup$
    – MarcoB
    Mar 9, 2017 at 6:03
  • $\begingroup$ @MarcoB I believe this belongs here, because - next to the theoretical aspect - it shows how to make use of Mathematica's abilities to simulate and to work with its distribution-framework. $\endgroup$
    – gwr
    Mar 9, 2017 at 12:34
  • $\begingroup$ My bad, I corrected my "decision support" (cf. my answer): In the simple 0-1-Utility case (right/wrong and a fixed amount to win), you would best choose $1$ (not 8, not 10). My exmple meets the case, when you will win an amount equal to the number of draws it took. $\endgroup$
    – gwr
    Mar 9, 2017 at 18:07

2 Answers 2

10
$\begingroup$

I would like to choose a slightly different approach here and point out, that this is an example for the Negative Hypergeometric Distribution. Wikipedia does give a nice table for this:

Choice of distributions

While for large numbers the distribution tends to the Negative Binomial Distribution the difference essentially is drawing with replacements (Binomial) vs. drawing without replacements (Hypergeometric).

Analytical Approach

In the present case, even when we may not have remembered all this statistical theory, we could still take an analytical approach by making use of ProbabilityDistribution. Let us look at the probabilities for the first few numbers of draws:

\begin{align} p(1) &= p(\text{ace}) =\frac{4}{52} \\ p(2) &= p(\neg\text{ace}) \times p(\text{ace}) = \frac{52-4}{52} \times \frac{4}{52-1} \\ p(3) &= p(\neg\text{ace}) \times p(\neg\text{ace}) \times p(\text{ace}) = \frac{52-4}{52} \times \frac{52-4-1}{52-1} \times \frac{4}{52-2} \\ \vdots \end{align}

So eventually we see a pattern here and can (hopefully) come up with the formula for the probability distribution function (in the discrete case a probability mass function) $f$ given the number of cards (or balls in an urn) $N$ and the number of marked cards (or balls) $M$:

\begin{align} f_{N,M}(k) = \frac{m}{N-k+1} \times \prod \limits_{i=1}^{k-1} \frac{N-M-i+1}{N-i+1} \end{align}

To implment this as a distribution we write the following code:

distribution[ n_Integer, m_Integer ] /; m <= n := ProbabilityDistribution[

    m/( n - k + 1) Product[ (n - m - i + 1)/(n - i + 1),{ i, 1, k-1 } ],
    { k, 1, n - m + 1, 1 }
    (* , Method -> "Normalization" is not needed here, but may be worth remembering
       for all cases where giving a proportional formula is easy to come up with *) 
]

analyticalDist = distribution[ 52, 4 ];

$PlotTheme = {"Detailed", "LargeLabels"};

Panel @ DiscretePlot[ 
    Evaluate@PDF[ analyticalDist, \[FormalX] ], 
    {\[FormalX], 0, 49, 1}, 
    ImageSize -> Large, 
    PlotLabel -> Style["Analytical PDF", "Subsection"]
]

Analytical PDF

Also we can find moments and what have you:

Mean @ analyticalDistribution (* e.g. the expected number of draws *)

$\frac{53}{5}$

Empirical Approach

We could of course also arrive at answers using simulation:

empiricalDist = EmpiricalDistribution @ Table[ 
    Min @ RandomSample[ Range[52], 4], 
    {10000} (* number of trials *)
];

N @ Mean @ empiricalDist

10.5809 (* true value again was 10.6 *)

If we compare the graphs for the PDF we note that the empirical distribution already is quite close:

Panel @ DiscretePlot[
    Evaluate @ {
        PDF[ empiricalDist, \[FormalX] ],
        PDF[ analyticalDist, \[FormalX] ]
    },
    { \[FormalX], 1, 49, 1 },
    PlotLegends -> {"empirical", "analytical"}, 
    ImageSize -> Large, 
    PlotLabel -> Style["PDF Comparison", "Subsection"]
]

Comparison of PDFs

General Case: Negative Hypergeometric Distribution

Grabbing a text book near you, you will find the formula for the probability distribution function of the above mentioned Negative Hypergeometric Distribution. It is a special case of the BetaBinomialDistribution - thanks to @ciao for pointing this out repeatedly until I listened :) - for those of you, who do not know this, there is a terrific paper showing the relations between univariate distributions by (Leemis and McQueston 2008) and more information can be found on Cross Validated.

Using this information, we can implement the Negative Hypergeometric Distribution explicitly as follows:

(* following Mathematica's parametrization for the hyergeometric case *)
negativeHypergeometricDistribution[k_, m_, n_ ] := TransformedDistribution[
    \[FormalX] + 1,
    \[FormalX] \[Distributed] BetaBinomialDistribution[ k, m, n-m ]
]

(* unfortunately parametrization is no standardized so one has to dwelve into
   the formulae for this in the paper, which are given though *)

nhg = negativeHypergeometricDistribution[ 1, 4, 52 ];

Mean @  nhg

$\frac{53}{5}$

And we can indeed convince ourselves that our case simply is a special case of a NHG-distributed random var $X$, where $X$ is the number of draws from $N = 52$ cards, where $M=4$ are "successes" and where we need $k=1$ success to stop:

And @@ Table[ PDF[ analyticalDist, i ] === PDF[ nhg, i ], {i, 1, 52 } ]

True

EDIT: Decision Support / Expected Utility

The OP has edited his question and asked for guidance in playing at a casino (e.g. what number of draws to bet on?). For this we have to turn to decision theory and excepted utility.

Let's assume that you are risk neutral. In the case that you get a fixed amount if you are right and zero if you are wrong, then it would be trivially optimal to bet on $1$. ($1$ has the greatest probability and the 0-1-loss function makes the expected utility equal to the probability; check this by simulation!)

Let us look at a more elaborate case to make the principle clear: If you were to receive the number of draws needed worth in money if you bet on the right number of draws, then your utility function would be $U(x,k) = x I(x = k)$ (where $I$ is the indicator function, $x$ your bet and $k$ the number of draws it took. In this case your best bet will be $13$:

NArgMax[ 
    {
        NExpectation[ 
          \[FormalX] * Boole[ \[FormalX] == i], 
          \[FormalX] \[Distributed] negativeHypergeometricDistribution[ 1, 4, 52 ]
        ],
        1<= k<= 49
    }, 
    k ∈ Integers
]

13

With[
    {
        dist =  negativeHypergeometricDistribution[ 1, 4, 52 ]
    },
    Panel @ DiscretePlot[ Evaluate@ Expectation[ 
          \[FormalX] * Boole[ \[FormalX] == i ], 
          \[FormalX] \[Distributed] dist 
        ], 
        {i,1,49}, 
        ImageSize -> Large, 
        PlotLabel -> Style["Expected Utility (0-1)", "Subsection"],
        PlotLegends->None
    ]
]

Expected Utility

$\endgroup$
4
  • 1
    $\begingroup$ "since it is not implemented...". It is, as BetaBinomial, already noted in an earlier reply. $\endgroup$
    – ciao
    Mar 9, 2017 at 23:35
  • $\begingroup$ @ciao Thanks, I did not make that connection. It is deeply burried in the docs and also not very obvious from web research. Guess one really needs to know what your looking for. :) $\endgroup$
    – gwr
    Mar 9, 2017 at 23:50
  • 3
    $\begingroup$ Not so much an MMA docs issue, more of a wider issue of distributions going by different names/definitions in the mathematical world. Would be nice if a future MMA docs edition (or MathWord) had a big interconnected chart with all alternate names/definitions and their relationships... $\endgroup$
    – ciao
    Mar 10, 2017 at 0:50
  • $\begingroup$ This is the best overview with regard to the relationships between univariate distributions I have found so far. $\endgroup$
    – gwr
    Mar 10, 2017 at 8:19
4
$\begingroup$

Update (attribution: @ciao)

As @ciao has pointed out, use ofBetaBinomialDistribution is the most straightforward for exact calculation and simulation. See the comment:

Mean[BetaBinomialDistribution[1, 4, 48]] + 1

or simulation

Mean[RandomVariate[Mean[BetaBinomialDistribution[1, 4, 48]] + 1],100000]
Median[RandomVariate[Mean[BetaBinomialDistribution[1, 4, 48]] + 1],100000]

etc

Original Answer

As the choice of first Ace is arbitrary you could simplify, e.g.:

tab = Table[Min[FirstPosition[RandomSample[Range[52]], #][[1]] & /@ Range[4]], 
   10000];
Histogram[tab, Automatic, "PDF"]
N@Mean[tab]
Median[tab]

enter image description here

You could also do analytically,e.g. let $f(k)$ be the probability that the first ace in the first $k$ cards, then:

f[k_] := 1 - Binomial[52 - k, 4]/Binomial[52, 4]
DiscretePlot[f[k], {k, 0, 52}]

Sanity checks:$f(k)=1$ for $k\ge 49$

The probability that is $k$th card is discrete derivative:

ListPlot[Differences[Table[f[k], {k, 0, 52}]], Filling -> Axis]
Differences[Table[f[k], {k, 0, 52}]].Range[52]

enter image description here

$\endgroup$
2
  • 1
    $\begingroup$ Or you could just Mean[BetaBinomialDistribution[1, 4, 48]] + 1... ;=} $\endgroup$
    – ciao
    Mar 9, 2017 at 6:56
  • $\begingroup$ @ciao doh! always forgetting to exploit diverse built-in distributions...hope your entrepreneurship is going well ;-) $\endgroup$
    – ubpdqn
    Mar 9, 2017 at 6:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.