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

  • $\begingroup$ 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. $\endgroup$
    – bzki
    Commented May 22, 2020 at 16:22
  • $\begingroup$ 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! $\endgroup$
    – lkjhgfdsa
    Commented May 24, 2020 at 8:35
  • $\begingroup$ 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$. $\endgroup$
    – bzki
    Commented May 27, 2020 at 13:46

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

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

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?

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

# 46.68259
# 2.89322

Somehow, the results are pretty normally distributed. enter image description here

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.