Counting Biased Coins 2: Estimating Bias Imbalance

Question

This is a continuation of another question I have asked recently. The setup and the question itself are sufficiently different (and likely more complicated) to warrant a separate question.

Setup: There is a device that produces coins. Most of the time it makes unbiased coins, but sometimes it makes biased coins. Let

• $p$ be the probability of producing a Heads-biased coin (P[Heads] > P[Tails])
• $q$ be the probability of producing a Tails-biased coin (P[Heads] < P[Tails])
• $1-p-q$ thus be the probability of producing an unbiased coin (P[Heads] = P[Tails])

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.

Questions:

1. We would like to test if $p > q$
2. We would like to estimate $p$ and $q$, and have some error bars or confidence regions.

Attempt: I have attempted a non-parametric approach.

1. For each coin, perform a binomial test in each direction
2. Count the number of significantly heads-biased and tails biased coins $N_H$ and $N_T$ given confidence threshold of $5\%$
3. Perform another binomial test on the two numbers, with the null hypothesis that the probabilities of producing a heads-biased and tails-biased coins are equal.

The good things about this solution is that it seems to look reasonable for actual data and that it gives us estimates for $\hat p = \frac{N_H}{N}$ and $\hat q = \frac{N_T}{N}$. The bad things are that I'm not not sure how solid the math is, the confidence threshold is somewhat arbitrary, and 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

In your last question, I provided a hierarchical Bayesian mixture model. That model can be generalized quite easily (to my surprise) to accommodate your changes here. The extension is not intended to be a legitimate solution to your problem, considering you're not familiar with Bayesian modelling. In any case, I'm posting it here for posterity.

I'll first present the full model and comment on its structure after. The model is

$$ \mathbf{p} \sim \operatorname{Dirichlet}(\mathbf{1}) $$

$$ \mu_l \sim \operatorname{Uniform(0, 0.5)} $$

$$ \mu_r \sim \operatorname{Uniform(0.5,1)} $$

$$ \kappa_l \sim \operatorname{Half Cauchy}(0,1)$$

$$ \kappa_r \sim \operatorname{Half Cauchy}(0,1)$$

$$ b_1 \vert \mu_r, \kappa_r \sim \operatorname{Beta}\Big(\mu_r \times \kappa_r, (1-\mu_r) \times \kappa_r \Big)$$

$$ b_2 \vert \mu_l, \kappa_l \sim \operatorname{Beta}\Big(\mu_l \times \kappa_l,(1-\mu_l) \times \kappa_l \Big)$$

$$ b_3 = 0.5 $$

$$ y_i \sim \sum_{i=1}^3 \mathbf{p}_i \operatorname{Binomial}(b_i;200) $$

The probability of drawing a biased coin ($p$ in the previous model, $\mathbf{p}$ in the present model) is now generalized to be a draw from a Dirichlet distribution. You can think of the elements of this vector as probabilities $p$, $q$, and $1-p-q$ from your description.

Each of the two biases are modeled as coming from a binomial distribution parameterized by the mean $\mu$ and the "precision" parameter $\kappa$. The means $\mu_l$ and $\mu_r$ are forced to be below and above 0.5 respectively.

The model is easy to write down, but challenging to fit in Stan. The Stan model is

data{
    int n;
    int y[n];
  }
  parameters{
    simplex[3] prob_of_bias;
    

    real<lower=0, upper=1> mu_left;
    real<lower=0> kappa_left;
    vector<lower = 0, upper = 0.5>[n] b_left; 
    
    real<lower=0, upper=1> mu_right;
    real<lower=0> kappa_right;
    vector<lower = 0.5, upper = 1>[n] b_right; 
    
    
    
  }
  model{
    real lp[3];
  
    prob_of_bias ~ dirichlet(rep_vector(1, 3));
  
    mu_left ~ beta(1, 1);
    mu_right ~ beta(1,1);
    
    kappa_left ~ cauchy(0, 1);
    kappa_right ~ cauchy(0, 1);
    
    // the prior below takes a parameterization in terms of mu and kappa
    // and turns it into alpha beta parameterization
    b_left ~ beta_proportion(mu_left, kappa_left);
    b_right ~ beta_proportion(mu_right, kappa_right);
    for (i in 1:n){
    lp[1] = log(prob_of_bias[1]) + binomial_lpmf(y[i] | 200, b_left[i]);
    lp[2] = log(prob_of_bias[2]) + binomial_lpmf(y[i] | 200, 0.5);
    lp[3] = log(prob_of_bias[3]) + binomial_lpmf(y[i] | 200, b_right[i]);
    target += log_sum_exp(lp);
    }
  }
   generated quantities{
  matrix[n,3] ps;
  for (i in 1:n) {
    vector[3] pn;
    // log-probability that there is bias
    pn[1] = log(prob_of_bias[1]) + binomial_lpmf(y[i] | 200, b_left[i]);
    // log-probability that there is no bias
    pn[2] = log(prob_of_bias[2]) + binomial_lpmf(y[i] | 200,  0.5);
    pn[3] = log(prob_of_bias[3]) + binomial_lpmf(y[i] | 200,  b_right[i]);
    // posterior probabilities for bias and no bias
    ps[i,] = to_row_vector(softmax(pn));
  }
   }

The model experiences divergences using sensible default parameters for Stan's samplers. I've had to decrease the step size of the integrator, thereby increasing the time required to sample the model. I can get the model to fit, but only just barely. My code could stand to be optimized. I make no attempt to optimize it at this time.

Simulation

I've simulated the process here. I have a 60% probability to drawing a fair coin, and a 25% probability to drawing a coin with bias less than 50%. Here is some code to generate some data:

# Simulate the data
ncoins = 100
nflips = 200
which_bias = sample(1:3, replace = T, size = ncoins, prob = c(0.25, 0.6, 0.15))
bias = matrix(rep(0, 3*ncoins), ncol = 3)
bias[,1] = rbeta(ncoins, 20, 80)
bias[,2] = 0.5
bias[,3] = rbeta(ncoins, 900, 100 )
theta = rep(0, ncoins)
for (i in 1:ncoins){
  theta[i] = bias[i, which_bias[i]]
}

y = rbinom(ncoins, nflips, theta)

Fitting the model with the simulated data results in the following marginal posterior distributions for the elements of $\mathbf{p}$. I've added the true values in red. Not bad for 100 coins and 200 flips.

Again, we can compute the posterior probability that any given coin is has a bias in a particular direction. Here are the probabilities first 10 coins

       P(Left Bias) P(No Bias) P(Right Bias)
    [1,]        0.000      1.000        0.000
    [2,]        0.000      0.000        1.000
    [3,]        0.001      0.999        0.000
    [4,]        0.000      1.000        0.000
    [5,]        0.002      0.998        0.000
    [6,]        0.000      1.000        0.000
    [7,]        0.000      0.000        1.000
    [8,]        0.000      1.000        0.000
    [9,]        1.000      0.000        0.000
   [10,]        0.000      1.000        0.000

To compare, the first 10 coins have the following biases:

 1 Fair      
 2 Right Bias
 3 Fair      
 4 Fair      
 5 Fair      
 6 Fair      
 7 Right Bias
 8 Fair      
 9 Left Bias 
10 Fair

So the probabilities are more are less bang on. That is only because $\mu_l$ an $\mu_r$ were very extreme biases.

Again, I'm not saying this is the solution you should use. It is a solution, and its a solution which could use a lot of optimization. If readers have extensions to the model or improvements, please feel free to add them.