Bernoulli trials

Sample size selection

1 Problem

We want to be able to select \(n\), the sample size, from a universe of Bernoulli trials with probability of success equal to \(\theta \in [0,1])\). For instance, consider a random sample \(x_1,\ldots,x_n\) from vodka bottles from a very large shipment intercepted by the police regarding potential contamination by methanol. Notationally, we would say that \[ x_1,\ldots,x_n \ \ \mbox{are iid} \ \ Bernoulli(\theta), \] or that \(Pr(x_i=1)=\theta\) and \(Pr(x_i=0)=1-\theta\).

2 Frequentist approach

As you probably have learned in an intro class on stats or probability, a \(100(1-\alpha)\%\) (approximate) confidence interval for \(\theta\) is given by \[ {\widehat \theta} \pm Z_{1-\alpha/2}\left[\theta(1-\theta)/n\right]^{1/2} \] where \({\widehat \theta}=(x_1+x_2+\cdots+x_n)/n\) and \(Z_{1-\alpha/2}\) is the \(100(1-\alpha/2)\%\) percentile of the standard normal distribution. If \(\alpha=0.05\), then \(Z_{1-\alpha/2}=1.959964\) (or simple the usual 1.96, or even 2!).

Therefore, the error has size \[ E = Z_{1-\alpha/2}\left[\theta(1-\theta)/n\right]^{1/2}. \] Since \(\theta(1-\theta)\) is always smaller than \(0.25\) (show it!), then \[ E \leq Z_{1-\alpha/2}0.5/\sqrt{n}. \] When \(\alpha=0.05\), then \(E \leq 1/\sqrt{n}\) and, consequently \(n \geq 1/E^2\). Suppose one wants the error to be at most \(E=0.01\), then \(n \geq 10,000\), while when \(E=0.05\), then \(n \geq 400\).

2.1 Selecting the sample size

For a given confidence level \(1-\alpha\), error size \(E\), it is easy to see that \[ n \geq \frac{Z_{1-\alpha/2}^2\theta(1-\theta)}{E^2}. \]

alpha  = 0.05
E      = 0.05
Z      = qnorm(1-alpha/2)
theta  = seq(0.01,0.2,length=100)
n      = Z^2*theta*(1-theta)/E^2

plot(theta,n,xlab="True proportion",ylab="Selected sample size")

3 Bayesian approach

Suppose that we do have some prior information about \(\theta\), represented here by a Beta distribution with parameters \(a\) and \(b\), ie. \(\theta \sim Beta(a,b)\). For instance, when \(b=15\) and \(a=b/9=1.67\), it is easy to see that \[ E(\theta)=0.1 \ \ \ \mbox{and} \ \ \ Pr(0.01 \leq \theta \leq 0.278) = 0.95. \] The 2008 Bayesian Analysis paper Bayesian Sample Size Determination for Binomial Proportions, by M’Lan, Joseph and Wolfson (Volume 3, Number 2, pp. 269–296, DOI:10.1214/08-BA310), offers various schemes to select the optimal sample size. From the paper’s abstract we collected the following part:

This paper presents several new results on Bayesian sample size determination for estimating binomial proportions, and provides a comprehensive comparative overview of the subject. We investigate the binomial sample size problem using generalized versions of the Average Length and Average Coverage Criteria, the Median Length and Median Coverage Criteria, as well as the Worst Outcome Criterion and its modified version.

For illustration, we use their result from Theorem 1, which uses the Average Length Criterion (ALC). to derive an approximate sample size formula, for when \(a>1\) and \(b>1\): \[ n = \frac{4Z_{1-\alpha/2}^2}{l^2}\left[\frac{B(a+1/2,b+1/2)}{B(a,b)}\right]^2-(a+b), \] where \(l\) is the length of the highest posterior density (HPD) interval.

alpha = 0.05
Z = qnorm(1-alpha/2)
b = 15
a = b/9
l = 2*E
n1 = 4*Z^2/(l^2)*(beta(a+1/2,b+1/2)/beta(a,b))^2-(a+b)

theta = seq(0,0.5,length=1000)
plot(theta,dbeta(theta,a,b),xlab="Proportion",type="l",ylab="Density")
abline(v=qbeta(0.025,a,b),lty=2)
abline(v=qbeta(0.975,a,b),lty=2)
abline(v=a/(a+b),lty=2)
title(paste("n=",round(n1),sep=""))

3.1 Varying the scale parameter \(b\)

Let us now see the behavior of \(n\) when \(b\), which is a scale parameter, varies.

alpha = 0.05
Z = qnorm(1-alpha/2)
b = seq(2,120,length=100)
a = b/9
l = 2*E
n11 = 4*Z^2/(l^2)*(beta(a+1/2,b+1/2)/beta(a,b))^2-(a+b)

par(mfrow=c(1,2))
plot(b,qbeta(0.025,a,b),type="l",ylim=c(0,1),ylab="95% prior credibility interval")
lines(b,qbeta(0.975,a,b))
abline(h=a/(a+b),lty=2)
abline(v=b[n11==max(n11)],lty=2)

plot(b,n11,type="l",ylim=c(0,100),ylab="Sample size")
abline(h=0,lty=2)
abline(v=b[n11==max(n11)],lty=2)

3.2 Highly informative prior

When your prior is highly informative, say \(\theta \sim Beta(120/9,120)\), whose prior mean is \(0.1\), the prior standard deviation is \(0.026\) and the prior 95% credibility interval is \((0.055,0.156)\), then you will only need a sample of size \(n=2\).

thetas = seq(0,0.5,length=1000)
plot(thetas,dbeta(thetas,120/9,120),type="l",xlab="Proportion",ylab="Prior density",lwd=2)
lines(thetas,dbeta(thetas,15/9,15),col=2,lwd=2)
legend("topright",legend=c("n=2","n=101"),bty="n",lwd=2,col=1:2)