Hypothesis test involving proportion of proportions

Question

I have 100 different coins, each with probability $p_i$ of being heads - but these $p_i$s are unknown.

I flip each coin $n_i$ times, and see that it lands on heads $x_i$ times. So the proportion of heads observed for each coin is $a_i := \frac{x_i}{n_i}$, so that $a_1, \dots, a_{100} \in [0,1]$.

Let $X$ be the number of $p_i$s that are $\geq 0.5$.

My null hypothesis $H_0$ is that $X=50$, i.e. half the $p_i$s are greater than or equal to $0.5$, and half are less than $0.5$.

Can anyone advise how I can test this hypothesis in R?

If needed, my prior assumption would be that $p_i=0.5$ for every $i$.

Answer 1

To set the ball rolling, we may sum up the cumulative posterior probability > 0.5 for each coin. For example, if one coin i flips ni=10 times with xi=3 heads, the chance of the true pi>0.5 would be the cumulative posterior probability > 0.5. In case of uniform prior, the posterior will have the exact shape of the Binomial likelihood yi = 1 - pbinom(0.5*10, size=10, prob=3/10) according to this answer.

The expected number of coins with pi>0.5, that is y=sum(yi), should be close to 50. A single simulation below returns y of 44.45 .

set.seed(1)
#single simulation
p <- c(runif(50, min=0, max=0.5), runif(50, min=0.5, max=1)) #exactly 50 and 50 with p<.5 and >.5
n <- sample(10, size=length(p), replace=TRUE) #number of draw for each ball
a <- rbinom(length(p), size=n, prob=p) / n #observed frequency

#adjust for extreme results
a[a==0] <- 0.05 / n[a==0]
a[a==1] <- 1 - 0.05 / n[a==1]

#expected number of p > 0.5
y <- function(a, n) {
  return(sum(pbinom(0.5*n, size=n, prob=a, lower.tail=F)))
}
y(a, n)
#44.44555

However, I'm stuck by the results of 100,000 simulations that show underestimated mean 46.68 . Maybe it's due to low number of coin flip xi<=10?

set.seed(1)
results <- NULL
for (i in 1:100000) {
  p <- c(runif(50, min=0, max=0.5), runif(50, min=0.5, max=1))
  n <- sample(10, size=length(p), replace=TRUE)
  a <- rbinom(length(p), size=n, prob=p) / n
  a[a==0] <- 0.05 / n[a==0]
  a[a==1] <- 1 - 0.05 / n[a==1]

  results = c(results, y(a, n))
}

mean(results)
# 46.68259
sd(results)
# 2.89322

Somehow, the results are pretty normally distributed.

# dev.new(height=4, width=4)
hist(results, breaks=100, probability=TRUE)
x_ <- seq(0, 100, by=0.1)
lines(x_, dnorm(x_, mean(results), sd(results)), col='red')