Counting Biased Coins

Question

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.

Answer 1

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

Answer 2

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 heads[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).

Your Questions

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.

Answer 3

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.