Problem 17 of Bernoulli's Ars Conjectandi

I'd like to solve Problem 17 of Bernoulli's Ars Conjectandi. I'll summarize the problem as:

A roulette wheel with 38 pockets labelled "1", ..., "8" (with 4 copies of each label), each pocket may hold at most 1 ball, and every pocket is equiprobable for a ball to land in, the player is given 4 balls. If the pocket labels are points, and the player sums the points awarded to them by the pockets the balls land in, what's the expected number of points the player may earn?

Solution

We can work out the frequencies we will find the balls producing k points.

For a single ball, it will land in a pocket labelled w with probability 4/32=1/8.

The second ball has two possibilities: it will land in a pocket also labelled w, or it will land in a pocket with a different label x. There are 31 empty pockets for the second ball to land in for both cases. But when the second ball lands in a w pocket, there are only 3 vacant w-labelled pockets (for a probability of 3/31). On the other hand, there are 4 vacant x-labelled pockets.

We can implement this in R as (similar to James Hanley's solution):

n <- rep(0,32)  # first 3 will remain 0, since points range is 4:32

for (label_1 in 1:8) {
  f1 <- 4; # 4 possibilites of a label_1
  for (label_2 in 1:8) {
    f2 <- f1 * (4 - (label_1 == label_2)); # 3 or 4 possibilities of label_2, depending...
    for (label_3 in 1:8) {
      f3 <- f2 * (4 - (label_3 == label_1) - (label_3 == label_2)); # etc
      for (label_4 in 1:8) {
        f4 <- f3 * (4 - (label_4 == label_1) -
                    (label_4 == label_2) - (label_4 == label_3));
        points <- label_1 + label_2 + label_3 + label_4
        n[points] <- n[points] + f4;
      }  
    }
  }
}

freq <- n[4:32]/24;

expected_value <- sum((4:32)*freq)/sum(freq); # = 18

In this algorithm, we need to divide through by 4! = 24 to avoid double counting. (Think of the case of getting 4 points, i.e., all four balls land in pockets with label "1": there's only one way for this to happen. We don't care how the result comes about [e.g., which ball landed in which particular pocket with label 1], we only care about the final configuration.)

Homework. Bernoulli gives us a "payoff table", rewarding the player with a number of coins depending on the points won. While computing the expected payoff amounts to sum(payoff*freq)/sum(freq), perhaps a better inverse problem is: determine what payoffs will result in the player's expected winnings to be exactly 4, such that if probability of points \(\Pr(p)\lt\Pr(p')\) the payoff \(f(p)\lt f(p')\) (or, \(f\) is monotonically decreasing on the interval [18, 32] and monotonically increasing on the interval [4, 18]).

(Historically, Bernoulli's payoff table was: 120, 100, 30, 24, 18, 10, 6, 6, 6, 5, 3, 3, 3, 2, 2, 3, 3, 3, 3, 4, 4, 6, 8, 12, 16, 24, 25, 32, 180. This is from 4 points up to 32 points.)

Lingering Puzzle

The problem people have with this worked example is, Bernoulli gives a solution for the expected winnings from his table as 4 + 349/3596 (approximately 4.09705228031). If we simulate the table with his winnings, we get a larger value (or smaller value, depending on our random number generator).

It's been a debate: has Bernoulli made a mistake? Professor E computed the expected winnings to be 4 + 153/17980 (about 4.00850945495), much lower than Bernoulli's answer.

An argument in defense of Bernoulli is, if we consider some code to generate the expected winnings:

simulate_wheel <- function(SIMS = 1000000) {
  Total <- 0
  nummi <- c(120, 100, 30, 24, 18, 10, 6, 6, 
             6, 5, 3, 3, 3, 2, 2, 3, 3, 3, 3, 
             4, 4, 6, 8, 12, 16, 24, 25, 32, 180)
  
  pocket_values = rep(1:8,4)
  
  INDICES <- 1:length(pocket_values)
  
  for(sim in 1:SIMS) {
    total_value <- sum(sample(pocket_values, 4, replace = F))
    Total <- Total + nummi[total_value-3]
  }
  
  Total/SIMS
}
simulate_wheel()

This produces a higher-than-expected amount of winnings consistently. (Your mileage may vary, depending on your random number generator.) But the expected sum of the labels for the pockets where the balls landed, we get:

simulate_wheel <- function(SIMS = 1000000) {
  Total <- 0
  
  pocket_values = rep(1:8,4)
  
  INDICES <- 1:length(pocket_values)
  
  for(sim in 1:SIMS) {
    total_value <- sum(sample(pocket_values, 4, replace = F))
    Total <- Total + total_value
  }
  
  Total/SIMS
}
simulate_wheel()

An expected 18 points, agreeing with earlier results.

Puzzle: Did Bernoulli make a mistake? Did our simulations mislead us? What's the real expected winnings according to the table given?

Homework 1. What's the expected winnings if a pocket could hold k balls (for k = 2, 3, 4)?