Comparing sequence of unfair coin flips with different predicted probabilities

Question

Suppose we have $n$ independent coins, each of which has an unknown and different probability of coming up heads. We have a magical machine that guesses the probabilities of each coin coming up heads. For $n=5$ coins, its guesses might look like 0.03, 0.45, 0.17, 0.91, 0.76. Then we flip each coin once and record the results, which might be heads, tails, tails, heads, heads. Our question: is there a statistical test that evaluates how accurate the machine is?

My first intuition was a chi-squared test like Pearson's or possibly a G-test, where our observed value is 1 if heads and 0 if otherwise. However, the data isn't categorical, nor do the observed values sum to a known constant. I also considered a Kolmogorov-Smirnov test to compare the expected with the observed, but neither are probability distributions.

I'm not sure how to phrase my preferences mathematically, but I'll try to do them qualitatively:

• edit: any prediction of, say, $p=0.5$ by the machine should be interpreted to say that the machine has evidence that 50% of coins with this probability will land heads and should be rewarded accordingly if the proportion of heads for coins with this probability is very close to 50%; it does not mean that the machine is clueless and ambivalent on heads/tails
• accuracy should be determined mostly on heads ratio: for instance, when the machine guesses 0.2 for every coin in a set, we should consider it accurate if 1 out of 5 landed heads or if 11 out of 50 landed on heads, but inaccurate if 3 out of 5 landed heads
• true probabilities might not be uniformly distributed from 0 to 1
• optional: accuracy should be reasonably tolerant of singular, highly unlikely events

Sorry if a similar question has been asked before—I'm very new to this site and statistics in general, but hopefully this question is clear enough.

Edit: As a response to a posted answer, we are not looking for Bernoulli proper scoring rules. To see why, suppose our machine predicts $p=0.5$ for all of 100 coins. In situation A, 50 coins land heads up, and in situation B, 95 coins land heads up. Clearly, the machine was a better predictor in situation A, which should be reflected in our test. However, the log-loss and Brier scoring rules, for example, each return the same score in both situations: $$\text{log-loss score, situation A}=-\frac{1}{100}\left(\sum_{n=0}^{50}\ln0.5+\sum_{n=51}^{100}\ln(1-0.5)\right)\approx0.69$$ $$\text{log-loss score, situation B}=-\frac{1}{100}\left(\sum_{n=0}^{95}\ln0.5+\sum_{n=96}^{100}\ln(1-0.5)\right)\approx0.69$$ $$\text{Brier score, situation A}=\frac{1}{100}\left(\sum_{n=0}^{50}(1-0.5)^2+\sum_{n=51}^{100}(0-0.5)^2\right)=0.25$$ $$\text{Brier score, situation B}=\frac{1}{100}\left(\sum_{n=0}^{95}(1-0.5)^2+\sum_{n=96}^{100}(0-0.5)^2\right)=0.25$$ The machine's predictions should be evaluated in aggregate. The machine's individual predictions should not be evaluated against its individual corresponding outcomes.

Answer 1

There are possibilities of formally testing the null hypothesis that your machine's predictions give the correctly calibrated probabilities of the experiments you observe.

Recall the definition of a p value. It starts with a test statistic, which is a single number derived from your observations. The p value is defined as the probability of observing a test statistics at least as extreme as the one you actually did observe, under the assumption that the null hypothesis is true.

What test statistic shall we use? The established way of assessing the quality of probabilistic predictions is proper scoring rules. I personally prefer the log loss, so that is what I will use here. An alternative would be the Brier score, see that last link.

The procedure is as follows:

1. Given your probabilistic predictions, generate a sequence of Bernoulli trials with the predicted success probabilities. For this generated sequence, calculate the log loss.
2. Do step 1 many times, e.g., 10,000 times. You now have 10,000 random log losses under the null hypothesis (because the observations were generated according to your predictions).
3. Calculate the log loss of your predictions for the actual observations you have.
4. Find out which proportion of the simulated log losses is larger than your actually observed log loss.

As an example, suppose that we have nine trials with actual success probabilities $q_1 = (0.1, \dots, 0.9)$. (Of course, we don't know this.) We have observed this sequence:

FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE  TRUE

First, assume that we actually do have predictions that match the true probabilities. Running the procedure above gives us the following histogram of log losses, with the actually observed log loss (of 0.591) indicated by a red vertical line:

We find that 2,977 of our 10,000 simulated log losses are larger (i.e., more extreme) than the actually observed log loss, so $p=0.2977$.

Second, assume that our prediction is wrong: instead of the true probabilities, we predict $q_2=(0.9, \dots, 0.1)$. In this case, we observe a log loss of 1.170. Under this null hypothesis, the log losses look like this:

This time, only 8 of 10,000 simulated log losses are larger than the actually observed one, so $p=0.0008$, and we will reject the null hypothesis that $q_2$ contains the true probabilities at standard alpha thresholds.

(If we had both predictions $q_1$ and $q_2$, then we could immediately say that $q_1$ is a better prediction than $q_2$ base on the lower log score alone. But that is not what you are asking.)

Of course, you can re-run the analysis using the Brier score or any other proper scoring rule for Bernoulli experiments.

Standard caveats of null hypothesis significance testing apply: statistical significance is not the same as clinical/economic/substantial significance; you can and will have both Type I and Type II errors; if your sample size is large, even tiny deviations from the true value will give you a significant result.

R code:

rm(list=ls())
probs_1 <- seq(0.1,0.9,by=0.1)
probs_2 <- 1-probs_1

set.seed(0)
(obs <- runif(length(probs_1))<probs_1)

log_score <- function(obs, preds) {
    -mean(obs*log(preds)+(1-obs)*log(1-preds))
}

sim_log_scores <- function(predictions, n_sims=1e4) {
    result <- rep(NA,n_sims)
    for ( ii in 1:n_sims ) {
        set.seed(ii)
        sim_obs <- runif(length(predictions))<predictions
        result[ii] <- log_score(sim_obs, predictions)
    }
    result
}

sims_1 <- sim_log_scores(probs_1)
(log_score_1 <- log_score(obs, probs_1))
hist(sims_1,main="Log scores under the null hypothesis probs_1",xlab="")
abline(v=log_score_1,lwd=2,col="red")
1-ecdf(sims_1)(log_score_1)

sims_2 <- sim_log_scores(probs_2)
(log_score_2 <- log_score(obs, probs_2))
hist(sims_2,main="Log scores under the null hypothesis probs_2",xlab="")
abline(v=log_score_2,lwd=2,col="red")
1-ecdf(sims_2)(log_score_2)