Three different approaches: Simulation, programming, and probability theory.
Riddles
Author

Jonatan Pallesen

Published

August 23, 2019

### Problem

A problem from the Riddler 5 Apr, 2019

Lucky you! You’ve won two gift cards, each loaded with 50 free drinks from your favorite coffee shop, Riddler Caffei-Nation. The cards look identical, and because you’re not one for record-keeping, you randomly pick one of the cards to pay with each time you get a drink. One day, the clerk tells you that he can’t accept the card you presented to him because it doesn’t have any drink credits left on it.

What is the probability that the other card still has free drinks on it? How many free drinks can you expect are still available?

### Simulation solution

Use Julia
library(JuliaCall)

I run a function use_until_empty where I subtract one from a random card until one of the cards is empty. At that point the function returns the number of drinks left on the other card.

A number of such simulations are run, and the results are saved in drinks_left

using Random
Random.seed!(0);
nsims = 10^6;
n = 50;

function use_until_empty()
cards = [n, n]
while minimum(cards) >= 0
cards[rand([1, 2])] -= 1
end
return maximum(cards)
end;

drinks_left = [use_until_empty() for _ in 1:nsims];

Probability other card has any drinks:

1 - sum(drinks_left .== 0) / nsims
0.92033

Number of free drinks you can expect are available (on average):

mean(drinks_left)
7.03872

### Dynamic programming

a is an array where position (i+1, j+1) is the probability that at some point you will have used i free drinks from one card and j free drinks from the other card.

ntickets = 50;
d = ntickets + 2;
a = fill(0.0, d, d);
a[2, 2] = 1;

for diag in 5:2*d, i in 2:d, j in 2:d
if i + j == diag
a[i, j] = a[i-1, j] / 2 + a[j-1, i] / 2
end
end;

d is the tile where 50 drinks are used from a card, so a[d, d] is the tile where 50 drinks have been used on both cards. The probability that the other card has any drinks left is 1 minus this value.

1 - a[d, d]
0.9204107626128212

The expected average number of free drinks available on the other card, can as above be found by summing over all probabilities a[i, d] (the chance that i free drinks have been taken on one card, while the other card is used) multiplied by 50 - i (the number of tickets that are left in those cases).

sum(a[i, d] * (d - i) for i in 2:d-1)
7.038512976105056

## Probability theory

Let the number of drinks on a card be n. Then let us consider what is the chance that the number of drinks taken from the other card is k.

We have to pick one card n times to empty it.

We have to pick the other card k times.

This is a total number of picks of n + k.

So it’s a from n + k choose n binomial.

$\binom {n + k} n (\frac{1}{2})^{n + k}$

One complication is that after this has been chosen, we also have to attempt to take another pick from the empty card, instead of the card that has been picked k times. The chance of this is 1/2. So we have to multiply the result by 1/2.

But also there are two ways to get a solution where one card has been picked n times and the other k times. Let us call one card A and the other B. One way is that card A is picked n times and then revealed to be empty. The other way is that card B is picked n times and then revealed to be empty. Since it can happen in two ways, we should multiply the probability by 2.

As it happens, this exactly cancels out with the first consideration, and we do not have to change the initial formula.

The probability that the other card has any drinks is 1 minus the chance of picking all 50 k:

n = BigInt(50);
k = BigInt(50);
1 - binomial(n + k, n) * (1/2)^(n + k)
0.9204107626128212385018729497578295385970684595752666786426521294828262398368679

The expected average number of free drinks available on the other card, can be found by summing over all probabilities k (the chance that k tickets have been taken) multiplied by 50 - k (the number of tickets that are left in those cases).

$\sum_{k=0}^{n} k \times \binom {n + k} n (\frac{1}{2})^{n + k}$

sum(binomial(n + k, n) * (1/2)^(n + k) * (50 - k) for k in BigInt(0):BigInt(50))
7.038512976105054911310832074459216601696085582898065457092134922234549776476342