Let $S \sim \text{Bin}(N, \pi)$ denote a number of successes. Using the non-informative $\text{Beta}(0.5, 0.5)$ prior, the posterior distribution of the probability of success is $$ \text{Beta}(0.5 + S, \, 0.5 + N - S) $$

Say that the interest is in the posterior probability that $\pi > 0.60$, i.e.

$$ \text{Pr}\big(\text{Beta}(0.5 + S, \, 0.5 + N - S) > 0.60\Big) $$

The R function f() below assumes that the true $\pi$ equals $0.8$ and computes the proportion of times, over 10000 simulations, that this posterior probability is $> 0.95$.

f <- function(N, pi=0.8)
  S <- rbinom(n=10000, size=N, prob=pi)
  proba <- 1 - pbeta(0.60, 0.5 + S, 0.5 + N - S)
  mean(proba > 0.95)

I would expect $f$ to be an increasing function of $N$ (after all, the posterior distribution becomes increasingly concentrated around $0.80 > 0.60$ as $N$ increases).


enter image description here

Do you have an explanation for the behaviour?


It occurs because $S$ is discrete.

For every value $N$, there is a certain value of $S$, call it $S_N$, such that the posterior probability of $\pi>0.6$ is greater than 0.95. We can find $S_N$ using this R function:

find_cutoff = function(N) {
  s = 0:N
  min(s[which(1-pbeta(.6, .5+s, .5+N-s) > .95)])

Now that we have the cutoff, we need to find the probability that $S\ge S_N|N$ which is just evaluating one minus the cdf of a binomial distribution.

d = ddply(data.frame(N=10:40), .(N), function(x){
  cutoff = find_cutoff(x$N)
  data.frame(cutoff = cutoff,
         prob = 1-pbinom(cutoff-1, x$N, .8))
plot(prob~N,d, type='b')

enter image description here


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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