# Counting Biased Coins

Edit Note: While this question is very interesting and relevant in its own right, I have come to a realisation that I have to make it a bit more complicated in order for it to be applicable to my actual research question. In order to avoid confusion, I have asked it as a follow-up question.

Problem:

There is a device that produces coins. Most of the time it makes unbiased coins, but sometimes it makes biased coins. Let $$q$$ be the probability of producing a biased coin. Biased coins can have different bias. There is no prior knowledge about the distribution of these biases. However, in practice we frequently see strongly-biased coins.

We use the device to produce $$N=100$$ coins. Then we toss each coin $$M=200$$ times and record the results. We would like to test if $$q > 0$$, with the null hypothesis that $$q=0$$. We would like to estimate $$\hat q$$ and get a confidence interval for it.

Clearly, the problem is unsolvable if we allow for coins of arbitrarily small but non-zero bias. In order to circumvent this problem, we will only be counting the strongly biased coins.

Given the bias of the coin $$b$$ we will define a coin as strongly-biased if $$|b-0.5|\geq 0.05$$. So the goal is to estimate the fraction of strongly biased coins $$\hat q'$$ that the device produces, as well as its confidence interval.

Attempt: I have attempted a non-parametric approach. It seems to work in practice, but I'm not sure if the mathematics behind it is solid, or if it is the most powerful test I can do.

1. For each coin, perform a binomial test against the null hypothesis, obtain a p-value
2. Threshold each of the obtained p-values to $$T = 5\%$$, obtain a list of 1s and 0s
3. Perform a binomial test on that list, with null hypothesis $$P[1]=T$$ and alternative hypothesis $$P[1] > T$$.

The logic behind it is that by testing multiple coins, some of them will test positive by chance. In order to answer the first question, we test if the number of coins we have classified as biased significantly exceeds that chance. My worry about this solution is that the threshold is somewhat arbitrary, that I don't know how it relates to the above definition of strongly-biased coins, and that I don't have a confidence interval.

It would be cool to alternatively try a model-based approach, addressing the hierarchical nature of this random process directly. However, I know nothing at all about the parametric methods used to tackle such problems. Any suggestions, names, links to literature are appreciated.

Note: This is a minimal example of a problem I have encountered in experimental design in neuroscience. I have judged that further details are unnecessary to make progress on this problem.

There are certainly many ways you could construct a test statistic giving the evidentiary ordering for this hypothesis test. Your method is coherent, but it is unlikely to be powerful. Since there is no a priori information on the degree of bias in a biased coin, a natural thing to do here would be to compute the maximum deviation from the expected value under fairness over all of the coins (or some simple scaled version of this). To do this, let $$X_{i,j}$$ denote the outcome of coin $$i$$ on toss $$j$$ and construct the test statistic:

$$R \equiv R(\mathbf{X}) \equiv 2 \times \max_{0 \leqslant i \leqslant N} \Bigg| \frac{1}{M} \sum_{j=1}^M X_{i,j} - \frac{1}{2} \Bigg|.$$

This test statistic is a scaled version of the maximum deviation from the expected value under fairness, where the scaling is used so that it takes on values in the range $$0 \leqslant R \leqslant 1$$. Larger values of the test statistic are more conducive to the alternative hypothesis that there is some non-zero probability of bias, so this gives a simple evidentiary ordering for the hypothesis test. Use of this evidentiary ordering gives the p-value function:

\begin{align} p(\mathbf{x}) &\equiv \mathbb{P}(R(\mathbf{X}) \geqslant R(\mathbf{x}) | H_0) \\[16pt] &= 1 - \mathbb{P}(R(\mathbf{X}) < R(\mathbf{x}) | H_0) \\[12pt] &= 1 - \prod_{i=1}^N \mathbb{P} \Bigg( M-R(\mathbf{x}) < \sum_{j=1}^M X_{i,j} < R(\mathbf{x}) \Bigg| H_0 \Bigg) \\[6pt] &= 1 - \prod_{i=1}^N \mathbb{P} \Bigg( M-R(\mathbf{x})+1 \leqslant \sum_{j=1}^M X_{i,j} \leqslant R(\mathbf{x})-1 \Bigg| H_0 \Bigg) \\[6pt] &= 1 - \Bigg( \sum_{s=M-R(\mathbf{x})+1}^{R(\mathbf{x})-1} \text{Bin}(s|M, \tfrac{1}{2}) \Bigg)^N. \\[6pt] \end{align}

(Note that the summation here is non-empty if and only if $$R(\mathbf{x}) \geqslant \frac{2}{M}$$, so if $$R(\mathbf{x}) < \frac{2}{M}$$ then we have a p-value of one.) Note that this test assumes that all coins are tossed the same number of times; if you were to remove this assumption then things would become more complicated and a generalised test statistic (accounting for the number of trials for each coin) would be desirable.

Programming the test: We can program this as a custom hypothesis test in R as follows. The input for the function is a binary matrix with rows corresponding to different coins (test objects) and columns corresponding to different trials.

binary.bias.test <- function(x) {

#Check validity of inputs
if (!is.matrix(x))             {
stop("Error: Input x should be a binary matrix") }
if (sum(x == 0) + sum(x == 1) != length(x)) {
stop("Error: Input x should be a binary matrix") }

#Set description of test and data
N           <- nrow(x)
M           <- ncol(x)
method      <- "Binary bias test"
data.name   <- paste0("data matrix ", deparse(substitute(x)), " with outcomes of ", N,
" test objects used for ", M, " trials each")
null.value  <- NULL
alternative <- "non-zero probability that test objects are \"unfair\""

#Calculate test statistics and p-value
RX          <- 2*max(abs(rowSums(x)/M - 0.5))
estimate    <- RX/2
statistic   <- RX
attr(estimate,  "names") <- "estimated maximum deviation from fairness"
attr(statistic, "names") <- "R"
if (statistic < 2/M) {
p.value <- 1    } else {
p.value <- 1 - (pbinom(RX-1, size = M, prob = 0.5) -
pbinom(M-RX, size = M, prob = 0.5))^N }

TEST        <- list(method = method, data.name = data.name,
null.value = null.value, alternative = alternative,
estimate = estimate, statistic = statistic, p.value = p.value)
class(TEST) <- "htest";
TEST }


Now let us examine this test using an example where we take $$q = 0.03$$ and we take $$p \sim \text{U}(0,1)$$ for a biased coin. As can be seen from the example, the test identifies the non-zero probability of bias in the coin, with a p-value of $$p < 2.2 \times 10^{-16}$$ for the test.

#Generate data
set.seed(1)
N <- 100
M <- 200
q <- 0.03
DATA <- matrix(0, nrow = N, ncol = M)
for (i in 1:N) {
BIASED <- (runif(1) <= q)
if (!BIASED) { DATA[i,] <- sample(c(0, 1), size = M, replace = TRUE) }
if (BIASED)  { PROB     <- runif(1)
DATA[i,] <- sample(c(0, 1), size = M, replace = TRUE,
prob = c(1-PROB, PROB)) } }

#Apply the test
binary.bias.test(DATA)

Binary bias test

data:  data matrix DATA with outcomes of 100 test objects used for 200 trials each
R = 0.77, p-value < 2.2e-16
alternative hypothesis: non-zero probability that test objects are "unfair"
sample estimates:
estimated maximum deviation from fairness
0.385

• This nested maximum is a bit difficult to decipher. I imagine that this might be more easy to write down. – Sextus Empiricus Jan 19 at 12:22
• Possibly --- I'm open to you proposing an edit (which I might reject!) if you have a preferable way of writing it. – Ben Jan 20 at 3:59
• I would write it by expressing it like the absolute value of the deviation (instead of the maximum of the negative/positive deviation) $$R \equiv R(\mathbf{X}) \equiv \max_{0 \leqslant i \leqslant N} \Bigg|\frac{1}{M} \sum_{j=1}^M (2 X_{i,j} -1)\Bigg|$$ I was also thinking about abbreviating the sum as a mean $$R \equiv R(\mathbf{X}) \equiv \max_{0 \leqslant i \leqslant N} \Bigg| (2 \bar{X}_{i} -1)\Bigg| \quad \text{with} \quad \bar{X}_{i} = \frac{1}{M} \sum_{j=1}^M X_{i,j}$$ but this notation for a mean with two indexes might not be so clear (are there alternatives?). – Sextus Empiricus Jan 20 at 9:54
• @SextusEmpiricus: Yes, that looks nicer --- edited. – Ben Jan 20 at 12:32
• I find this maximum a smart way to tackle the problem. The distribution of the maximum absolute deviation will shift to larger values in both cases that the size of the bias increases or that the frequency of the biased coins increases. But, I wonder whether possibly a different order statistic might be more powerful. For instance, if the frequency is high but the bias is small then possibly using the second largest deviation might be more sensitive. Possibly some scheme like the Holm Bonferroni method could be more powerful. – Sextus Empiricus Jan 20 at 14:53

There are a few solutions from a Fequentist perspective, but I'd like to offer a Bayesian perspective.

## The Model

It would be useful to model the data generating process as a mixture of binomials; one with 0 bias and the other with some non-zero bias left to be estimated. You can then partially pool the data from all the coins to help estimation. Here is the model I'm using in math...

$$\mu \sim \operatorname{Beta}(1,1)$$ $$\kappa \sim \operatorname{Half-Cauchy}(0,1)$$ $$p \sim \operatorname{Beta}(1,1)$$ $$b_i \sim \operatorname{Beta}(\mu \times \kappa, (1-\mu) \times \kappa)$$ $$y_i \sim p \operatorname{Binom}(n_i, b_i) + (1-p)\operatorname{Binom}(n_i, 0.5)$$

Here, $$p$$ is the probability of producing a biased coin, $$b_i$$ is the probability of heads for the biased coins, which is assumed to come from a beta distrbution with parameters $$\alpha = \mu \times \kappa$$ and $$\beta= (1-\mu) \times \kappa$$.

If I were to write this in Stan, it might look like

data{
int n;
int coin[n];
}
parameters{
real<lower=0, upper=1> prob_of_bias;
vector<lower = 0, upper = 1>[n] b;

real<lower=0, upper=1> mu;
real<lower=0> kappa;
}
model{
prob_of_bias ~ beta(1,1);
mu ~ beta(1,1);
kappa ~ cauchy(0, 1);
// the prior below takes a parameterization in terms of mu and kappa
// and turns it into alpha beta parameterization
b ~ beta_proportion(mu,kappa);
for (i in 1:n){
target+=log_mix(prob_of_bias,
binomial_lpmf(heads[i] | 200, b[i]),
binomial_lpmf(heads[i] | 200,  0.5));
}


Here is a simulation example. The probability of generating a biased coin is 0.77. Here is the simulation code

n = 100 #Number coins
is_biased = rbinom(n, 1 , 0.77)
b = 0.5 + is_biased*(rbeta(n, 90, 10) - 0.5)
x = rbinom(n, 100, p) #100 flips
d = tibble(coin = factor(1:n), heads = x)


Fitting this model, we see the estimated posterior means of each parameter recover their true values to within 1 decimal place (which is fine because the credible intervals capture the truth).

The question about the fraction of strongly biased coins has now been turned into a question of estimating $$p$$ from the model -- which is the probability of drawing a biased coin. We can obtain a variety of credible intervals but perhaps the easiest would be to just compute quantiles of the posterior for $$p$$.

Additionally, we can examine the probability that any given coin is biased by computing the posterior probability directly. Appending the following generated quantities block to the Stan model

 generated quantities{
matrix[n,2] ps;
for (i in 1:n) {
vector[2] pn;
// log-probability that there is bias
pn[1] = log(prob_of_bias) + binomial_lpmf(heads[i] | 200, b[i]);
// log-probability that there is no bias
pn[2] = log1m(prob_of_bias) + binomial_lpmf(heads[i] | 200,  0.5);
// posterior probabilities for bias and no bias
ps[i,] = to_row_vector(softmax(pn));
}


Will produce a matrix in which the rows correspond to coins and the columns are probabilities of being biased and not being biased respectively.

As an example, the first few coins from my simulation are

52 83 90 85 59


The middle 3 are biased coins. The first 5 rows of the matrix ps rounded to 4 decimal places are

      P(bias) P(No bias)
[1,]       0          1
[2,]       1          0
[3,]       1          0
[4,]       1          0
[5,]       0          1


In this case, the probabilities are quite large, but when the bias is smaller these probabilities might be more reasonable (e.g. 30% and 70%). In such a case, you have a probabilistic model for bias and can leverage techniques from that literature (or from Bayesian decision theory, setting some sort of utility function) to decide if a coin is biased or not.

I don't expect you to actually use this approach, and I'm sure similar approaches can be made with Expectation Maximization. Just thought I would throw my in two cents in case someone else found it interesting. Additionally, if this is something you are interested in, I can add details as needed.

• Thanks a lot, this is very interesting. I will need a bit of time to digest what is going on. I see you use some kind of hierarchical model, where parameters themselves are random variables. That sounds highly relevant, but I don't know at all how to work with such constructions. Could you hint me to the name of this construction (is it simply "hierarchical bayesian model"?), and perhaps lecture notes or any place where I can see a basic example solved with as much detail as possible? Thanks again for your effort – Aleksejs Fomins Jan 20 at 23:47
• Yes, "hierarchical bayesian model" is the name of the construction. Chapter 5 of Bayesian Data Analysis of Gelman et. al is a good resource. – Demetri Pananos Jan 20 at 23:49
• Ok, I understand now that if I follow this route, I will eventually want to arrive at the marginal posterior $P[p|Data]$, where $p$ is the fraction of biased coins. This will involve computing some complicated integrals over the joint distribution. I don't know Stan at all. To what extent can it simplify my life for this task. I could write all expressions explicitly in python and use some Gibbs sampler to do the integrals numerically. Is this procedure automatized in Stan? How hard is the installation / how steep is the learning curve? – Aleksejs Fomins Jan 21 at 11:17
• @AleksejsFomins You've understood correctly. I'm not sure how much you could simplify the task. Gibbs samplers take a long time to converge to the posterior, so learning Stan is your best bet. Alternatively, pymc3 is a library in python which does the same thing. Installation is easy (see cmdstanpy library), but the learning curve is steep for the first little bit. – Demetri Pananos Jan 21 at 15:47
• @AleksejsFomins I can try to help if you'd like – Demetri Pananos Jan 21 at 15:47

I am not sure just what your question is asking you to do, Here are some relevant computations from R:

Test $$H_0:p=.5$$ vs $$H_a: p < .5$$ at level 5% (or lower).

If $$p = .5,$$ then $$P(X \le 41) = 0.0443,$$ where $$X$$ is the number of Heads in 100 tosses. So you can have a one-sided test at level 4.4%. Changing the critical value to 42 would give a test above the 5% level.

pbinom(41, 100, .5)
[1] 0.04431304
pbinom(42, 100, .5)
[1] 0.06660531


If $$p = .45$$ then $$P(X \le 41)= 0.2415.$$ So the power of the test at level 4.4% is only about 24% Thus, only about 24% of coins biased with $$p = 0.45$$ will be detected as biased (using critical value 41).

pbinom(41, 100, .45)
[1] 0.2414906


By contrast, if the critical value 41 is used, a coin biased with $$p=0.3$$ will be detected a little more than 99% of the time.

pbinom(41, 100, .3)
[1] 0.9928264


Thus, for a test at level just below 5%, 100 tosses is enough to do a pretty good job of detecting coins with bias $$p = .45,$$ but not so good at detecting biased with \$p = .45.%

Note: If you use normal approximations for the binomial probabilities, you might delude yourself into thinking you have a test at exactly the 5% level, but that is an artifact of the continuous normal approximation:

$$\mu = np = 50,\; \sigma = \sqrt{np(1-p)} = 5.$$

qnorm(.05, 50, 5)
[1] 41.77573


For the discrete binomial distribution, "critical value" 41.78 amounts to integer 41.