# Hypothesis test involving proportion of proportions

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$$.

• A little confused on what's the assumption in the problem, and what is being tested? Do we assume that $p_i = p$ for all $i$, and are we then setting as our null hypothesis that $p = 0.5$ (it would follow that $E(X) = 50$)? Alternatively, do we assume that $p_i = p$ but we don't assume $p = 0.5$, and our null is that $E(X) = 50$? It might help to clarify some of these items. Interesting problem.
– bzki
May 22, 2020 at 16:22
• Rereading my question it was very unclear - I was confusing my $p_i$s and $a_i$s. $X$ is actually the number of $p_i$s that are $\geq 0.5$, while I had previously written $a_i$s. I've rewritten the question - let me know if this is any better! May 24, 2020 at 8:35
• I have to think this over a bit. I will try to come back to this at some point this week if no one else has come up with something. My initial thoughts, which could be wrong, is that $X$ (despite the notation suggesting a random variable) is some fixed parameter that is the sum of fixed parameters $I(p_i \geq 0.5)$. So we would probably test the null hypothesis that $X = 50$ based on a random variable, call it $Y$, that is the sum of the observed proportions being greater than 0.5, i.e., $I(a_i \geq 0.5)$. I guess the trick is to know the null distribution of this under $X = 50$.
– bzki
May 27, 2020 at 13:46

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

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')