Geometric Mean for Combining Probabilities, addendum

We've noted before the geometric mean is best when combining several different forecasts together into one. Today I'd like to discuss how to do this specifically for election forecasts.

If we naively try to combine the forecasts that, e.g., Biden will win Nevada, using the forecasts as of August 30th:

Forecaster	Pr(Biden)	Pr(Trump)
DecisionDeskHQ.com	77.8%	22.2%
The Economist	89%	11%
FiveThirtyEight	77%	23%
JHKForecasts.com	86.1%	13.9%
OurProgress.org	79%	21%
PluralVote.com	76.6%	24.4%
ReedForecasts.com	62.5%	37.5%

We would obtain \[ \Pr(Biden) = 77.87021\%\tag{1a} \] but we would also find \[ \Pr(Trump) = 20.33664\%\tag{1b} \] ...which sums to 98.20685%, which is odd. These probabilities should always sum to 100%, what gives?

The solution is to first transform the probabilities into odds. Then take the geometric mean of the odds, and finally transform back to probabilities.

Why does this work? Well, as odds, a forecast would look like \[ O(\text{Biden wins NV}) = \frac{\Pr(Biden)}{\Pr(Trump)} = \frac{N_{\text{Biden}}}{N_{\text{Trump}}}. \tag{2}\] Taking the geometric mean of the odds gives us better approximations for the ratio of frequencies, which could then be transformed back into probabilities. We can obtain the probability of Biden winning Nevada from the odds by \[ \Pr(Biden) = \frac{O(Biden)}{1 + O(Biden)} \tag{3}\] and for Trump we could note \[ O(Trump) = 1/O(Biden) \] then find \[ \Pr(Trump) = \frac{O(Trump)}{1 + O(Trump)} = \frac{1}{1 + O(Biden)}. \tag{4}\] Hence we find adding Eq (3) to Eq (4) that \[ \Pr(Biden) + \Pr(Trump) = 1 \tag{5} \] probabilities sum to 100%, as expected and desired.

Applied to our forecasts, we find the odds given as:

Forecaster	Odds(Biden)
DecisionDeskHQ.com	3.504505
The Economist	8.090909
FiveThirtyEight	3.347826
JHKForecasts.com	6.194245
OurProgress.org	3.761905
PluralVote.com	3.273504
ReedForecasts.com	1.666667

The geometric mean of these odds is approximately 3.829059, hence a probability of Biden winning Nevada approximately 79.29203% and Trump has a 20.70797% chance of winning Nevada. This makes a difference for Biden of about 2%, whilst negligible improvement for Trump.

Puzzle/Homework. Consider the case of, say, a primary with several candidates. Suppose we have multiple forecasters make predictions for each candidate to win the primary. How can we generalize our method to handle this case?

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)?

Fast Dirichlet Distribution Parameter Estimation

From the Bayesian perspective, looking at \(\operatorname{Dir}(\alpha_{1}, \dots, \alpha_{n})\) the parameters \(\alpha_{j}\) encode the amount of evidence or number of observations of rolling an n-sided die, we witness side j approximately \(\alpha_{j}\) times.

It is tempting to use the Dirichlet distribution for estimating voting behavior, because a voter casts their ballot for a Democrat, a Republican, or a third-party candidate. The problem with this is, if we just set \(\alpha_{1}\) to the number of votes for the Democratic candidate (and likewise for the other parameters), there is too little variation. If we've seen a die turn up 1 and 6 roughly 99% of the time, we can be confident it will do so in the future. (We have such confidence about the Sun rising, after all.)

The problem: given observed outcomes \(x_{j}\) (dice rolls), we can form the approximate probabilities \(p_{j} = x_{j}/\sum_{\ell}x_{\ell}\). Using these probabilities, we can estimate the parameters of the Dirichlet distribution. I'm really interested in finding sensible values for \[\alpha_{0} = \alpha_{1} + \dots + \alpha_{n}.\tag{1}\] This quantity in Eq (1) is called the Precision of the Dirichlet distribution.

More precisely, we will observe an election in K counties, so we will end up with a vector \(\vec{x}_{j}\) whose components are the votes received by candidate j in each county. This gives us a vector \(\vec{p}_{j}\) whose components are the proportion of the county's vote went to candidate j

Solution: we can make estimates of the Dirichlet distribution based on the data, then use this to estimate the various parameters of interest.

Different Algorithms

There are a variety of ways to estimate parameters, I'm looking at closed-form algorithms. In particular, the method of moments derived algorithms.

In a nutshell, the method of moments looks at the moments of random variables \(\mathbb{E}[X^{n}]\) in terms of the parameters of the distribution, then tries to solve for one of the parameters in terms of the expected values. When given a population, we substitute in the sample mean (and sample variance) for the expected values in these formulas, giving us numerical estimates.

Minka's Method

Minka (2012) gives the following scheme. The method of moments applied to the Dirichlet distribution with parameters \(\alpha_{i}\), if we let \(\alpha_{0} = \alpha_{1}+\dots+\alpha_{n}\), may be approximated from data \(\vec{p}_{i}\) (in the trials the proportion corresponding to \(\alpha_{i}\)) using \[\widehat{\alpha_{0}} = \frac{\mathbb{E}[p_{1}] - \mathbb{E}[p_{1}^{2}]}{\mathbb{E}[p_{1}^{2}] - \mathbb{E}[p_{1}]^{2}}\] we find the approximation \[\boxed{\widehat{\alpha_{i}} = \mathbb{E}[p_{i}]\widehat{\alpha_{0}} = \mathbb{E}[p_{i}]\frac{\mathbb{E}[p_{1}] - \mathbb{E}[p_{1}^{2}]}{\mathbb{E}[p_{1}^{2}] - \mathbb{E}[p_{1}]^{2}}}\tag{2}\] However, there is no reason why we should use \(p_{1}\) in determining \(\widehat{\alpha_{0}}\) instead of any other \(p_{i}\).

Note algebraically speaking: \[\mathbb{E}[X^{2}] = \operatorname{var}[X] + \mathbb{E}[X]^{2}\]

A simple R implementation:

minka_estimate <- function(data) {
    mu <- mean(data[,1]);
    s2 <- var(data[,1]);

    (mu*(1 - mu))/s2 - 1;
}

Ronning's Method

We can use all the data to estimate \(\alpha_{0}\), as Ronning (1989) suggested. We use \[\operatorname{var}[p_{j}] = \frac{\mathbb{E}[p_{j}](1 - \mathbb{E}[p_{j}])}{1 + \alpha_{0}}\] after some algebraic rearrangement, then taking the log of both sides \[\log(\alpha_{0}) = \log\left(\frac{\mathbb{E}[p_{j}](1 - \mathbb{E}[p_{j}])}{\operatorname{var}[p_{j}]}-1\right).\] We take the average of the right-hand side over j (the candidates) to get \[\boxed{\log(\alpha_{0}) = \frac{1}{n}\sum^{n}_{j=1}\log\left(\frac{\mathbb{E}[p_{j}](1 - \mathbb{E}[p_{j}])}{\operatorname{var}[p_{j}]}-1\right).}\tag{3}\] The astute reader will recognize this is the geometric mean of estimates.

As a sanity test that this makes sense, we recover the correct value for the Beta Distribution.

A simple R implementation:

ronning_estimate <- function(data) {
    assert_that(all(0 <= data) && all(data <= 1));
    assert_that(all(0.97 <= rowSum(data)));
    
    col_estimate <- function(col) {
        mu <- mean(col);

        if (mu == 0) return(0);
        
        log(mu*(1-mu)/var(col) - 1);
    }
    
    exp(mean(apply(data, 2, col_estimate)));
}

Crazy algebraic results

With 3 results (a "trinomial", if you will), we have an estimate I found by hand. Let \(\mu_{A}\), \(\mu_{B}\) and \(\mu_{C}\) be the sample means and \(\sigma^{2}_{B}\) the sample variance. Assuming \[ 0\lt \mu_{A}\lt 1,\quad\mbox{and}\quad 0\lt\mu_{B}\lt 1-\mu_{A}\] (and if this fails to hold, permute the labels "A", "B", "C" until it holds, if at all possible), then we can estimate \[\widehat{\alpha_{A}} = \left(\frac{(1 - \mu_{B})\mu_{B}}{\sigma^{2}_{B}} - 1\right)\mu_{A} \tag{4.1}\] \[\widehat{\alpha_{B}} = \widehat{\alpha_{A}}\frac{\mu_{B}}{\mu_{A}} \] \[\widehat{\alpha_{C}} = \widehat{\alpha_{B}}\frac{1 - \mu_{B}}{\mu_{B}} - \widehat{\alpha_{A}}\]

Observe this uses Minka's formula in estimating \(\widehat{\alpha_{0}}\) (the parenthetic term on the right hand side of Eq (4.1)) for use in computing \(\widehat{\alpha_{A}}\). The reader can verify if we add these quantities together, we in fact get Minka's estimate using the second largest sample mean (and its associated sample variance).

Testing Results

When testing which estimator to use, we should have some error implicit in our mind, in the sense of "What are we trying to optimize? What counts as 'wrong'?" It's one thing to try to recover the parameters of a Dirichlet distribution, it's another to estimate the concentration parameter for a heuristic.

We could note that the method of moments for the Beta distribution coincide with estimating \(\alpha_{0}\) for using Minka's method on any particular outcome (i.e., using \(\widehat{\alpha_{j}}\) for any particular, fixed j). It's not uncommon for Bayesians to use a Gamma distribution as a prior for the concentration parameter of a Beta distribution, which would give us some measure of its variance.

The method of moments is biased (slightly). Despite this, we can still use the Cramer-Rao bound to estimate the variance of our estimates.

For what its worth, when used on voting results, the method of moments produce similar estimates when using the proportion of votes the Republican and Democrat received. It's within the margin of error of the maximum likelihood estimates given by MASS::fitdistr, but is consistently smaller than maximum likelihood estimates. Since I'm advocating humility in making these estimates, I prefer the method of moments using the fancy-pants algebraic method I derived (Minka's method on the runner-up candidate's proportion of votes).

If I were less lazy, I'd include more detailed proofs about estimates and goodness-of-fit tests, but...meh, exercise for reader :)

References

1. Jonathan Huang, "Maximum Likelihood Estimation of Dirichlet Distribution Parameters". Unpublished manuscript(?), undated, Eprint
2. Thomas P Minka, "Estimating a Dirichlet distribution". Microsoft Research Paper, 2012. Eprint.
3. G. Ronning, "Maximum-likelihood estimation of dirichlet distributions". Journal of Statistical Computation and Simulation 32 (1989) 215–221.

Puzzles from Bernoulli

I recently stumbled across a fascinating puzzle from Bernoulli's Ars Conjectandi, and then found that part 3 of his book consists of 24 equally exciting worked problems. But I haven't found online the presentation of these problems.

See either Anders Hald's A History of Probability and Statistics and Their Applications before 1750 (2005), especially chapter 15, section 5...or Edith Dudley Sylla's translation The Art of Conjecturing, Together with Letter to a Friend on Sets in Court Tennis (2006).

If I misrepresented any of the problems, leave me a comment! I was rather quick and rushed in assembling the exercises, and I easily could have made mistakes. Also, I am fully aware some of these problems are ambiguous. Try working out every variant you can think of. The whole point (at least, to me) is that these present variations on a theme.

Problem 1. There are two balls in an urn, a "winner" ball and a "loser" ball. There are three players. The first player draws a ball. If it's the "loser" ball, the first player returns it to the urn; and if it's the "winner" ball, the first player wins the game. Should the first player lose, the second player performs the same task, and wins only if drawing the "winner" ball. Should the first and second players both draw the "loser" ball, the third player draws a ball. If the third player draws the "loser" ball, the house wins the game. What is the expected winnings for each player (and the bank)?

Problem 2. A variant of problem 1, each player bets some amount. If no player draws the winning ball, then they divide up the bets equally. (Example: player 1 bets 4 coins, player 2 bets 2 coins, player 3 bets 1 coin; if they all draw the "loser" ball, then each player receives (4 + 2 + 1)/3 coins back.) What is the expected winnings for each player?

Problem 3. Consider a tournament (i.e., a sequence of 2-player games). Two players compete in a game, where there will be only one winner. The winner plays the game against player 3. The winner of the second match plays against player 4. And so on until player 6. If the lots of each player in the games are equal, how do the lots of the later players compare/increase over those of the earlier players?

Problem 4. As a variant of problem 3, the lot of the player who wins the first game is stipulated to be double that of the third player, who plays only in the second game, and so on. Compare the expected winnings for each player.

To be clear, Bernoulli means the relative probability mass of winning the first round are even for both players (each player has 1 favorable outcome for an odds of 1:1, or 50% probability of winning), but the relative probability mass for the winner of the first round to beat player 3 is doubled (the victor of the first match has twice the favorable outcomes as compared to game 1 [i.e., 2 favorable outcomes] whereas player 3 has 1 favorable producing a 2/3 probability favoring the victor of the first game to win the second), and the relative probability mass of the winner of the second round is doubled to determine the odds of winning the third game (so if the victor of the first game also won the second game, this victor has 4 favorable outcomes, but player 4 has 1 favorable outcome, the probability of winning this particular game for the victor of the first two games would be 4/5), and so on.

Problem 5. "This is the third problem of Huygens' Appendix." A wagers against B, that of 40 cards, of which 10 of each color, he will draw 4 of them in a way to have one of each color. (Or, for modern readers, consider a standard playing deck with Jacks, Queens, Kings discarded; A wages that B cannot draw one card from each suit.) What is the probability that A wins the bet? [Alleged solution: One finds in this case that the chance of A is to that of B as 1000 is to 8139.]

Problem 6. "This is the fourth problem of Huygens' Appendix." One takes 12 tokens of which 4 white and 8 black. A wagers against B that among 7 tokens that he will draw from them blindly, there will be found 3 white. One demands the ratio of the chance of A to that of B; i.e., what are the odds A wins? (There is some ambiguity in the historic text; namely, is it exactly 3 white or at least 3 white?)

Face Cards

Problem 7. Let there be a single face card in a pack of n cards, and the first player (of m players) to draw it wins. If no one draws it, and cards remain, then the players continue drawing cards. What is the probability of each player to win? What if n = mk (the number of players divides the number of cards)?

Problem 8. As a variant of Problem 7, what if there were j face cards in the deck and the winner is still the first person to draw a face card? What is the probability for each player to win?

Problem 9. As a variant of Problem 8, what if the players keep drawing cards until the deck is exhausted, and the winner is the player with the most face cards drawn? For ties, the winnings are split among the winners. (Assume each player pays 1 coin to play, for example.) What is the expected winnings for each player?

Problem 10. As an extension to Problem 9, what if we permit any player to sell their position? Specifically Bernoulli considers four players (A, B, C, D) with a deck of 36 cards, of which 16 are face cards. Each player receives cards in rounds until 23 cards have been distributed, A has received 4 face cards, B has 3, C received 2, and D a modest 1 face card, so that there remain 13 cards among which there are 6 face cards. The fourth player D (who is next to receive a card), "seeing that almost all hope of his winning has vanished", wishes to sell his right to one of the others. How much should he sell it for and what are the expectations of the individual players?

Dice Puzzles

Problem 11. Throw a die. Then throw a second die. If they differ, the player wins a point; if they agree, the player loses a point. Then the player throws a third die. If it agrees with any previous die, lose one additional point for each agreement; otherwise, the player wins an additional point added to their running score. Do this for a total of 6 dice. What is the expected score for the player?

Problem 12. Similar to Problem 11, but the dice must be thrown in numerical order. E.g., the first die must be "1" for the player to get a point, otherwise the player loses a point; the second die must be "2" for the player to get a point, otherwise the player loses a point; and so on. What is the expected score for the player?

Problem 13. Three players (A, B, C) have a list of 6 numerals ("1", "2", ..., "6") on a sheet of paper before them. Each of them take turns round-robin. On a player's turn, they roll a die, and if the result is on their sheet of paper, then the player gets to eliminate it from his sheet (scratch off the number) and roll again; but if the player has already eliminated that number, then the next player gets the die. This process continues until someone eliminates all their written numerals. It happens, however, after a while that A has 2 numerals before him, B has 4, and C has 3; it's A's turn to throw. What are probabilities for each player to win? [Bernoulli notes This problem requires more labor and patience than ingenuity.]

Problem 14. There are k players. A given player throws a die, which shows its face to be m. This tells the player to throw m dice and sum the values shown. (We may choose m to be added to the score or not.) The player with the most points wins.

Here's a twist, though: one player may opt to beat a fixed number t of points. This value t is fixed by the rules of the game.

What's the expected winnings for each player if no one opts for the fixed points route? What's the expected winnings for each player if someone chooses to beat a fixed number of points?

Although Bernoulli didn't pitch it, what if there's a "bidding war" competition to determine which player may opt to beat a fixed number of points? Instead of having t be fixed, when one player asks to be the one to beat a fixed number of points, that player must offer a bid ("I want to beat x points"). Each player may pass (and no longer participates in the bidding anymore), or offers a higher bid ("I want to beat y > x points"). This continues until the maximum value of 35 is bid. What strategy works best in this bidding strategy?

Problem 15. As a variant of Problem 14, what if the fixed number of points is the square of the first toss of the die?

Problem 16 (Cinq et Neuf). This is a prototype of craps. Player 1 tosses a pair of dice. Player 1 wins if on his first toss she gets a 3, an 11, or any pair. But player 2 wins if player 1 tosses a 5 or 9.

If player 1's first toss is a 4, 6, 7, 8, or 10, then the game continues (player 1 keeps tossing a pair of dice) until either (a) a 5 or 9 appears [player 2 wins], or (b) player 1 tosses the first value [player 1 wins].

What is the probability of player 1 winning? (Player 2's probability of winning is, by definition, the complement of the probability for player 1 winning.)

Although not asked, what is the expected number of tosses for a given game?

Wheel Games

Problem 17. The player pays 4 coins to toss 4 balls on a roulette-like wheel. The roulette wheel has 32 pockets, each with labels "1", "2", ..., "8". There are four pockets labelled "1"; four pockets labelled "2"; and so on. Each pocket may contain at most 1 ball. The player wins the sum of the pockets's labels (for the pockets containing the balls). What is the expected winnings for the player? [Solution.]

Card Games

Problem 18 (Trijaques). We consider a "toy model" of poker.

We assemble a deck of 28 cards from a standard playing deck, by discarding the cards 2 through 8 for each suit. The values for each card are determined by its face value, but the Jack of Clubs and 9s are wild.

The player will be given 4 cards. The goal is for the player to assemble either a flush (a run of all 4 cards regardless of suit, e.g., "9, 10, J, Q"), or a pair, three of a kind, or four of a kind. The value of a hand is the sum of the value of the cards in the combination. The player with the highest valued hand wins the pot (or it is split among the highest valued hands).

But the sequence of play is as follows: each player is dealt 2 cards face down. Then the players bet. Then each player is dealt 2 cards face up.

What is the expected winnings for each player? What strategies could be considered in the betting process?

Problem 19. Consider a generic game, where one player is the "banker". The banker has an advantage over the other players (i.e., is more probable than any other player to win a given round). But the rules of the game may allow moving the banker role to another player.

Specifically, the banker has probability p of winning a round, and probability q of losing, with p + q = 1 and \(r = p - q > 0\). The banker has probability h of continuing the next round as banker, and probability k of losing the position as banker to another player, with \(t = h - k > 0\). Let a denote the amount won by either the player or the banker, whoever wins the round.

What is the expected winnings for the banker after "many" rounds?

Problem 20 (Capriludium, Bockspiel). At the beginning of the game, each player puts down their bet. Then the banker shuffles the deck, and divides it into equally sized hands. Each player (and the banker) gets a hand. Bernoulli says the player just turns the hand over without organizing it. If the punter [player who is not the banker] has their facing card be of equal or higher value compared to the banker, then the punter wins an amount of money equal to his bet/bid from the banker. Otherwise the player loses their wager to the banker. When the banker loses to all the players in a game, the next player becomes the banker.

After one round, the top cards are not yet discarded. New wagers are first made. Then the top card is discarded (collected by the banker for later shuffling).

If there are N = sf cards in the deck with s suits and f face cards (of value ranging between 1 to f), and suppose there are n players (including the banker) for n = 2, 3, 4.

What is the expected winnings for the banker? What is the expected number of rounds to be played for a given deck? How does it vary on the number of players? What is the probability the banker will remain in their role as banker after one hand? After h hands?

Problem 21 (Basset). The basic formalization is there are 2n cards, of which k are marked "a" and \(2n - k\) marked "b". For example, 2n = 52, and k = 4 (e.g., aces in a standard deck of playing cards). The player draws two cards in succession (no replacement). The possibly outcomes:

• ab = the banker wins 1 point
• ba = the player wins 1 point
• aa = the banker wins 1 point
• bb = toss the cards aside and draw 2 new cards (and consult this table of outcomes again)

What is the expected winnings for the banker after 1 hand? For 1 game (exhausting the whole deck)?

As a variant, we could try using the full rules for Basset, with the bizarre bet multiplying schemes.

Curious Puzzles

Problem 22. There are two players, Titus and Caius. Titus pays Caius one coin for each round where Titus will throw a single die. There are a possible outcomes, of which b favor Titus (Caius pays him one coin) and c favor Caius (where Titus wins nothing). If Titus throws one of the c cases continuously n times in a row, Caius must return all n coins to Titus. What is the expected winnings of Caius and Titus?

Problem 23 (Blinde Würffel, "Blind Dice"). We have 6 dice with a number on only one face, and blanks on the remaining faces. One die has "1" for its non-blank side, another has "2" for its non-blank side, and so on, so each label "1" through "6" may possibly show up. Blank sides are treated as having value 0. Suppose the player rolls all 6 dice, and wins the sum of the values shown. What is the player's expected winnings?

Problem 24. A variant of Problem 23, if the player gets no points, in 5 tosses in a row then the player gets his money back for those 5 tosses. What is the player's expected winnings now?

2016 Exit Polls

Exit polls, by their nature, are extremely noisy and do not provide margins of error. Until 2016, news organizations banded together to provide a single, coherent exit poll for presidential elections. This coalition started collapsing after the 2016 presidential election. Both CNN and Fox News have slightly different exit poll data for the 2016 election. How do we make sense of these polls? When can we say they tell us "the same story"?

Toy Problem: Coin Tossing

Lets consider a simpler toy problem: I flip a coin N₁ times and obtain y₁ heads ("successes") and you flip a coin N₂ times and obtain y₂ heads ("successes"). How do we know our coins are equally biased? Assuming each of us has done a "large number of trials" (a couple hundred each).

The null hypothesis: these coins follow the same distribution with probability p of success.

The alternative hypothesis: these coins do not follow the same distribution.

Let \(p_{1} = y_{1}/N_{1}\) (and similarly \(p_{2} = y_{2}/N_{2}\)), a better approximation to the probability of heads ("success") is \[p=\frac{y_{1}+y_{2}}{N_{1}+N_{2}}.\tag{1}\] We can standardize the data and approximate it as a normal distribution with standard deviation \(\sigma_{1}^{2} = p(1-p)/N_{1}\) and similarly \(\sigma_{2}^{2} = p(1-p)/N_{2}\). Then the "true variance" is \[\sigma^{2} = \sigma_{1}^{2} + \sigma_{2}^{2} = p(1-p)\left(\frac{1}{N_{1}}+\frac{1}{N_{2}}\right)\tag{2}\] which we should use for performing the Z-transform: \[Z_{1} = \frac{p_{1} - p}{\sigma}\tag{3}\] and a similar definition of Z₂. The test statistic we want to measure is \[z = Z_{1}-Z_{2} = \frac{p_{1}-p_{2}}{\sigma}\tag{4}\] which should follow a normal distribution with mean 0 and standard deviation 1.

We can then pick some confidence level, usually 1 - α = 95% confidence level, which leads to the critical value \[z_{1 - \alpha/2} = \Phi^{-1}(1 - \alpha/2)\tag{5}\] using the quantile for the normal distribution \(\Phi^{-1}\) for approximately 1.96 for the 95% confidence level.

If the quantity computed in Eq (4) is "more extreme" than the quantity computed in Eq (5), i.e., if \(z_{1-\alpha/2} < |z|\), then we reject the null hypothesis and conclude these coins follow different distributions. Otherwise, we fail to reject the null hypothesis. Note: we never prove they follow the same distribution, we just fail to prove they follow different distributions.

Exit Polls

We can apply this same setup to exit polls. Unfortunately for Fox News, contrary to appearances, they do not provide state-level exit poll data for the first 12 questions. Perhaps this is a technical error due to some software bug on their server-side, but I can't do anything to remedy the problem.

Further, the results are identical percentages (on the first dozen questions) for FOX and CNN, despite having apparently different sample sizes. It is not hard to prove (if these results are accurately reported) the results of computing Eq (4) for each category is 0 identically, and thus we fail to reject the null hypothesis for each question.

It's rather anti-climactic despite this long buildup, but it's good to know these exit polls are coherent, in some appropriate sense.

The exit polls are in CSV form on github.

What is...a Random Variable?

Tentative Definition. A random variable assigns to a given "random phenomenon" or "random process" some (real) number, or a vector ("list") of numbers.

Examples. The following are all random variables.

1. Let X be the number of heads in 10 coin flips.
2. Let R be the number of times a given baseball pitcher strikes out an opponent in the course of a given game.
3. Let S be the number of "successes" (heads) until the first "failure" (tail) in a repeated trial (flipping a coin over and over again).
4. Let T be the waiting time (in minutes) until the next bus arrives.

Slightly more formal Definition. If we represent the possible outcomes for a given "random process" (i.e., the set underlying its sigma algebra) by Ω, then a random variable is a function \(X\colon\Omega\to\mathbb{R}\) (or possible \(X\colon\Omega\to\mathbb{R}^{k}\) some [fixed positive integer] k) such that for any \(x\in\mathbb{R}\) the preimage of smaller values is an event \(\{\omega\in\Omega:X(\omega)\leq x\}\in\Sigma\) ("is measurable", i.e., we can assign a probability to that preimage). We will write \(X\leq x = \{\omega\in\Omega:X(\omega)\leq x\}\) in an abuse of notation.

Remark. We should remember a sigma algebra is not just the space of outcomes \(\Omega\), but also the specific set of well-defined events \(\Sigma\subseteq\mathcal{P}(\Omega)\). There are various specifications we have on the events, for instance "something happens" \(\Omega\in\Sigma\); for any event \(E\in\Sigma\) its complement is also a well-defined event \(\Omega\setminus E\in\Sigma\) [or "an event does not happen"]; for any countable family of events \(\{E_{j}\}_{j\in J}\subseteq\Sigma\) their union is also a well-defined event \(\bigcup_{j\in J}E_{j} \in\Sigma\) ["one of these events might happen"]; and so forth. Implicitly we consider a probability measure on the set of well-defined events generically denoted \(\Pr(-)\). Altogether, this is the data necessary to describe some "random phenomenon".

Example. Let E be any event. Then indicator function \(I_{E}\) which is zero for any \(x\notin E\) and \(I_{E}(e)=1\) for all \(e\in E\), this is a random variable. For multiple events \(E_{1},\dots, E_{n}\), we find \(I_{E_{1}\cup\dots\cup E_{n}}(x)=\max(I_{E_{1}}(x),\dots,I_{E_{n}}(x))\) and \(I_{E_{1}\cap\dots\cap E_{n}}(x) = \prod_{j} I_{E_{j}}(x)\). For discrete random phenomena, these indicator functions are the basic building blocks for constructing other random variables.

What happens to all that sigma-algebra baggage? Given a random variable, we can ask "What is the probability its value will be in a given range?" For example, "What is the probability the starting pitcher for the Dodgers will strike out at least 30 batters?"

This would be computed by first assembling all possible outcomes which satisfy this \(\mathcal{E}=\{E\in\Omega : R(E)\geq30\}\subseteq\Omega\) and using the probability measure on the sample space for the random process \(\Pr(\mathcal{E})\) or more imaginatively \(\Pr(R\geq30)\).

As a caveat, we hasten to add, only inequalities are necessarily well-defined, e.g., \(\Pr(X\leq x)\). Equality "on the nose" may not be well-defined \(\Pr(X=x)\), but we abuse notation to write the Probability mass function in this manner. (This gets tricky with subtle nuances when dealing with continuous random variables instead of discrete ones.)

This induces a nice mathematical structure on the image of the random variable \(R(\Omega)\), namely we can "transport" the probability distribution from the sigma algebra \(\Sigma\) on \(\Omega\) to \(R(\Omega)\). This is the Distribution Function for the random variable, \(F_{R}(x) = \Pr(R\leq x)\).

Equivalence relation of Random Variables. If we have two random variables, say, \(X\) and \(Y\), we can say they are Equivalent if for any \(x\in\mathbb{R}\), we have \(\Pr(X\leq x) = \Pr(Y\leq x)\). This is usually denoted \(X\sim Y\).

Algebra of Random Variables. Given some random variables X, Y on the sample sigma algebra, we can define new random variables \(X+Y\), \(X - Y\), \(XY\), \(X/Y\) provided Y is never zero, and exponentiation \(X^{Y} = \exp(Y\log(X))\). The intuition to have is that the operations are done as operations on real-valued functions.

So, specifically, if \(x\in\Omega\), then \((X+Y)(x)=X(x)+Y(x)\), \((X - Y)(x) = X(x) - Y(x)\), \((XY)(x) = X(x)Y(x)\), \((X/Y)(x) = X(x)/Y(x)\) provided \(Y(x)\neq0\), and exponentiation \((X^{Y})(x) = \exp(Y(x)\log(X(x)))\).

Probability Distributions. We have a few "standard" distributions which are the "template" for various random processes. Flipping a coin follows a Bernoulli Distribution if we flip the coin only once, and a Binomial Distribution if we flip it N times, for example.

The notation for these families may vary reference to reference. A Bernoulli distribution with probability p of success is usually denoted \(\mathrm{Bernoulli}(p)\) or \(\mathrm{Ber}(p)\).

To indicate a random variable is distributed like one of these standard distributions, we abuse notation and write \(X\sim \mathrm{Bernoulli}(p)\).

We can build more distributions out of a handful of basic ones, for example \(Y = X_{1} + \dots + X_{N}\) where the \(X_{j}\sim\mathrm{Bernoulli}(p)\) will describe flipping a coin N times and counting the number of "successes" ("heads"). This gives us the Binomial distribution when we consider \(\Pr(Y\leq k)\) (there are at most k heads in N coin flips).

We can specify a probability distribution by its parameters (like p in the Bernoulli distribution), the probability mass function and/or the probability density function. Often it's useful to give other summary statistics alongside this data.

Expected Value. We also have for any random variable X its Expected Value given by \(\mathbb{E}[X] = \sum_{x\in X(\Omega)}x\Pr(X=x)\) (or replacing the sum with an integral for continuous random variables). The intuition we should have for the expected value of a random variable is this captures the "average value" of the random variable.

If we have some function \(f\colon\mathbb{R}\to\mathbb{R}\), then we have \(\mathbb{E}[f(X)] = \sum_{x\in X(\Omega)}f(x)\Pr(X=x)\) (and again, an integral instead of a sum for continuous random variables, with the restriction that f is an integrable function).

Note that \(\mathbb{E}[X^{2}]\neq(\mathbb{E}[X])^{2}\) and more generally \(\mathbb{E}[XY]\neq \mathbb{E}[X]\mathbb{E}[Y]\). But we do have \(\mathbb{E}[X + Y] = \mathbb{E}[X] + \mathbb{E}[Y]\) and, for any real number \(a\in\mathbb{R}\) \(\mathbb{E}[aX] = a\mathbb{E}[X]\).

Exercise. Let \(E\in\Sigma\) be an event in a sigma algebra, and \(I_{E}\) be the indicator function on E. What is \(\mathbb{E}[I_{E}]\)?

Theorem. Let X and Y be random variables, a and b be real numbers. Then:

1. \(\mathbb{E}[aX + b] = a\mathbb{E}[X] + b\)
2. \(\mathbb{E}[X+Y] = \mathbb{E}[X] + \mathbb{E}[Y]\)
3. \(\displaystyle\mathbb{E}[XY] = \sum_{\omega\in\Omega}X(\omega)Y(\omega)\Pr(\omega)\)
4. \(\displaystyle\mathbb{E}[X/Y] = \sum_{\omega\in\Omega}\frac{X(\omega)}{Y(\omega)}\Pr(\omega)\) if  \(Y(\omega)\neq0\) for any \(\omega\in\Omega\)
5. \(\displaystyle\mathbb{E}[X^{Y}] = \sum_{\omega\in\Omega}\exp(Y(\omega)\log(X(\omega)))\Pr(\omega)\)

Variance. If expected value tells us what neighborhood a random variable is likely to live in, the variance tells us how spread out this neighborhood is. We define it as \[\mathrm{Var}[X] = \mathbb{E}[X^2] - (\mathbb{E}[X])^2 = \sum_{\omega}(X(\omega) - \mathbb{E}[X])^{2}\Pr(\omega).\] The variance and expected value for a random variable contain a lot of useful information, which we use when trying to infer parameters from data.

More generally, we have for any two random variables X and Y a measure of their failure to be correlated by the covariance \[\mathrm{Cov}(X,Y) = \sum_{\omega}(X(\omega)-\mathbb{E}[X])(Y(\omega)-\mathbb{E}[Y])\Pr(\omega)\] which is such that \(\mathrm{Cov}(X,X)=\mathrm{Var}[X]\). Correlatedness is not the same as independence: independent random variables are not correlated, but uncorrelated random variables may or may not be independent (a "all salmon are fish, not all fish are salmon" type statement). So uncorrelated is a "weaker" property than independence.

Applications?

Iterate! Note, we can transform the parameters of these distributions (like p the probability of success in a Bernoulli trial) into random variables themselves. This is precisely what Bayesian data analysis does: the parameters are random variables following prior probability distributions, which we update as new data becomes available using Bayes's theorem.

Regressions! A linear regression basically says that the observations are really values of a random variable, i.e., \(Y\sim\mathcal{N}(aX + b, \varepsilon)\) where \(\mathcal{N}\) is the normal distribution. There are other useful regressions, but this is the basic idea.

We should admit that this is one formulation of regressions in terms of random variables. The other uses conditional random variables, \(Y|X\sim f(\beta\cdot X,\theta)\) when the regressions (X) are stochastic (i.e., "not controlled by the experimenter/statistician"). Formally these are different models. But when actually doing the regressions, they are treated "the same".

Tests! We often do an experiment, producing some data points \(x_{1},\dots,x_{n}\) which we interpret as values of a random variable X which follows a prescribed distribution. We test the assumption (that X follows the given distribution with specific parameters) by comparing the sample mean \[ \mu = \frac{1}{n}\sum_{j=1}^{n}x_{j} \] to the expected value \(\mathbb{E}[X]\). The central limit theorem suggests that \(\sqrt{n}(\mu - \mathbb{E}[X])\) looks like a normal distribution centered at 0 with variance approximately equal to the variance of the data points (loosely speaking).

Reading List

• Discrete Random Variables and their Expected Values, MIT OCW 18-05 Introduction to Probability and Statistics
• Expected Value and Variance, lecture notes.
• Statistical Machine Learning, by Han Liu and Larry Wasserman, specifically chapter 12 on Bayesian inference using random variables.

How many news stories are there?

Recently, an eccentric billionaire bought the Los Angeles Times and sought to make it rival the New York Times as a "newspaper of record". Presumably this means hiring more journalists, but let us ask a simpler question.

Puzzle 1: How many news stories go unreported by both the New York Times and the Los Angeles Times?

We can solve this puzzle using the maximum likelihood estimator for the Hypergeometric Distribution. Think of it like this: on a remote island with some unknown deer population, we go and (without harming the wildlife) tag K deer. A month later, we return, and capture n deer, of which k are tagged. We can estimate the total population of deer N on the island.

Explicitly connecting that analogous problem to our own, we know the "tagged stories" K reported by the New York Times, the "sample stories" n reported by the Los Angeles Times, of which there is the "tagged sample stories" k reported by both newspapers, and we want to estimate how many news stories there are in total N. The maximum likelihood estimator for N is given by \[ \min_{\widehat{N}}\frac{\Pr(\widehat{N},K,n,k)}{\Pr(\widehat{N}-1,K,n,k)}\geq1 \] the smallest N for which the ratio of probabilities is greater than 1. It is not hard to solve this to find \(\widehat{N} = [Kn/k]\) where the brackets indicate we are using the integer part of the number (e.g., [3.2]=3, [4.9]=4).

Now we just need to list the stories which the New York Times reported but the Los Angeles Times did not (giving \(K-k\)), the stories which both papers reported (k), and the total number of stories the Los Angeles Times reported (n). From this, we will estimate how many stories have gone unreported.

To answer this fully, I looked at the front section for each paper for May 28, 2019. The short answer is K = 22 stories in the New York Times, n = 12 stories in the Los Angeles Times, and k = 5 stories in both. We thus may expect there to be N = [264/5] = 52 stories, of which 29 were reported and 23 went unreported by either newspaper. Find below a density plot of the probability for various N, and notice how it is maximized at N = 52 (indicated by a red vertical line):

Solution: Using the maximum likelihood estimate for the hypergeometric distribution, there were a total of N = 52 news stories, 29 were reported by one of the two newspapers, and 23 stories went unreported.

Puzzle 2: Is there a Bayesian estimate for the number of news stories? Or different ways to estimate the total number of news stories?

Puzzle 3: How stable is this estimate for N? If we examine, say, the last week's worth of articles, do we get approximately the same value for N?

Puzzle 4: What if we extend this analysis to include, e.g., the Wall Street Journal, the Washington Post, and others? How stable is N in this case?

Find two tables below, one listing the stories in the international section for both papers, and the second for national stories. Corresponding stories are listed on the same row.

New York Times	Los Angeles Times
She Thought She’d Married a Rich Chinese Farmer. She Hadn’t. (A4)
Attacks by Extremists on Afghan Schools Triple, Report Says (A4)
Romania’s Most Powerful Man Is Sent to Prison for Corruption (A6)
With Trump’s Visit to Japan, Empress Masako Finds a Spotlight (A8)
Trump and Abe’s ‘Unshakable Bond’ Shows Some Cracks in Tokyo (A8)	Trump pushes off war talk on Iran, says ‘regime change’ is not U.S. goal (A1)
Election Puts Europe on the Front Line of the Battle With Populism (A10)	In European vote, far-right surge fails to materialize, but mainstream parties lose support (A2)
European Parliament Elections: 5 Biggest Takeaways (A10)
European Vote Reveals an Ever More Divided France (A11)
18 Schoolchildren Stabbed, and Girl and Man Killed, in Attack in Japan (A11)	Knife-wielding man attacks schoolgirls in Japan, killing 2 (blurb of story on A2)
Sebastian Kurz, Austrian Leader, Is Ousted in No-Confidence Vote (A12)	Ousted by parliament, Austria’s Kurz vows to win back job (A4)
Israel’s Netanyahu Struggles to Form a Government, as Time Runs Short (A12)	Netanyahu running out of time to form government; Israel may face new elections (A2)
White Panda Is Spotted in China for the First Time (A12)
30 Dead and 200 Missing in Congo After Boat Sinks (A12)
	Arrests, killings strike fear in Thailand’s dissidents: ‘The hunting has been accelerated’ (A3)

Matches are based on substantially overlapping subject matters. The only debatable story match is "Trump pushes off war talk on Iran", which is a proper subset of the corresponding New York Times article.

Also note, in the Los Angeles Times, there was a 1000 word blurb about the knife attacks in Japan. Later, on their website, they posted a longer and more detailed article. I decided to count that as a match, which may be debatable.

Sources: Los Angeles Times, New York Times

The national stories in both newspapers, appears to be completely disjoint sets of stories.

New York Times	Los Angeles Times
Trump Administration Hardens Its Attack on Climate Science (A1)
Google’s Shadow Work Force: Temps Who Outnumber Full-Time Employees (A1)
Trump Wants to Wall Off Huawei, but the Digital World Bridles at Barriers (A1)
With His Job Gone, an Autoworker Wonders, ‘What Am I as a Man?’ (A1)
With the 2020 Democratic Field Set, Candidates Begin the Races Within the Race (A1)
Saving Charlie: A Rush to Rescue Stranded Cats and Dogs from Oklahoma Floods (A17)
Fearing Supreme Court Loss, New York Tries to Make Gun Case Vanish (A17)
A Missed Opportunity for the Malpractice System to Improve Health Care (A19)
Why a Hamptons Highway Is a Battleground Over Native American Rights (A22)
	High radiation levels found in giant clams of Marshall Islands near U.S. nuclear dump (A1)
	He made millions as an L.A. investor. Now, he may run for president to fight poverty (A1)
	Want to park in Koreatown? Get ready for a ‘blood sport’ (A1)
	Put your hands together for the World Series of Poker, turning 50 this year (A4)
	Texas lawmakers approve safe gun storage program, quietly going around the NRA (A4)
	Oklahoma’s opiod lawsuit targeting drugmaker goes to trial Tuesday (A7)

Matches are based on substantially overlapping subject matters.

Sources: Los Angeles Times, New York Times

Geometric Mean in Probability

If we have N estimates for the probability of an event, say p₁, ..., p_N, then the best a good estimate for the probability is the geometric mean:

p

=

[p1×...×pN]1/N.

To see this, think of probability from a frequentist perspective, p_j = (estimated number of trials where event occurred)/(estimated number of trials), i.e.,

pj

=

nj,x / nj

where n_j,x is the number of trials where the event x occurred, and n_j is the estimated number of trials.

We can best estimate the numerator as the geometric mean of the numerators of our estimates

nx

=

[n1,x×...×nN,x]1/N.

(and similarly for the denominator), since the geometric mean is the best way to combine different estimates of possibly different orders of magnitude.

If the numbers are "close enough", the geometric mean will not differ greatly from the arithmetic mean ("average"). To see this, simply consider the Taylor expansion of [1 + x]^1/N to linear order, take p_j to be μ + Δp_j where μ is the arithmetic mean and |Δp_j /μ| < 1. Expanding the geometric mean will produce the arithmetic mean plus "small" corrections of order N⁻¹.

For numbers which are "spread far apart", the geometric mean gives better estimates than the arithmetic mean. I suppose one way to think about this is the logarithm tells us the "order of magnitude" for a quantity. The order of magnitude for the revised estimate should be the arithmetic mean of the orders of magnitudes for our various estimates. The geometric mean, as the revised estimate, is the only quantity that can do this.

One fun book (among many) is Order of Magnitude Physics.

Addendum (April 22, 2019). I struck out "the best" estimator, because I actually don't have a proof off the top of my head that this is optimal in any sense. It's "good", consistent with the frequentist interpretation of probability, and most importantly it works. But I do not have a well-defined notion of an "error" or "loss function" which the geometric mean of probability estimates minimize, and thus I felt it dishonest to describe it as "the best".

That said, "absence of evidence is not evidence of absence". I may be ignorant of some folklore that the geometric mean of probabilities optimizes some desirable property, and really is (in some sense) "the best estimator". I just don't have the proof to back the claim, so I will revise the claim.