a motivating example

Suppose I have two trick coins:

• one (coin A) comes up heads 75% of the time, and
• the other (coin B) only 25% of the time.

Possible answers:

1. Er, probably coin (A)?

2. Well, $$\begin{array}{rl}\mathbb{P}\{\text{6 H in 10 flips}|\text{coin A}\}& =(\genfrac{}{}{0ex}{}{10}{6})(.75{)}^6(.25{)}^4\\ & =0.146\end{array}$$ and $$\begin{array}{rl}\mathbb{P}\{\text{6 H in 10 flips}|\text{coin B}\}& =(\genfrac{}{}{0ex}{}{10}{6})(.25{)}^6(.75{)}^4\\ & =0.016\end{array}$$ … so, probably coin (A)?

For a precise answer…

3. Before flipping, each coin seems equally likely. Then

   $$\begin{array}{rl}\mathbb{P}\{\text{coin A}|\text{6 H in 10 flips}\}& =\frac{\frac{1}{2}\times 0.146}{\frac{1}{2}\times 0.146+\frac{1}{2}\times 0.016}\\ & =0.9\end{array}$$

Bayes’ rule

$$\begin{array}{r}\mathbb{P}\{B|A\}=\frac{\mathbb{P}\{B\}\mathbb{P}\{A|B\}}{\mathbb{P}\{A\}},\end{array}$$

where

• $B$: possible model
• $A$: data
• $\mathbb{P}\{B\}$: prior weight on model $B$
• $\mathbb{P}\{A|B\}$: likelihood of data under $B$
• $\mathbb{P}\{B\}\mathbb{P}\{A|B\}$: posterior weight on $B$
• $\mathbb{P}\{A\}$: total sum of posterior weights

Suppose instead I had 9 coins, with probabilities 10%, 20%, …, 90%; as before I flipped one 10 times and got 6 heads. For each $\theta$ in $0.1,0.2,\dots ,0.8,0.9,$ find $$\begin{array}{r}\mathbb{P}\{\text{ coin has prob }\theta |\text{ 6 H in 10 flips }\}.\end{array}$$

Question: which coin(s) is it, and how sure are we? (And, what does it mean when we say how sure we are?)

Uniform prior

Weak prior

Strong prior

The likelihood: 6 H in 10 flips

The likelihood: 30 H in 50 flips

The likelihood: 60 H in 100 flips

The likelihood: 6,000 H in 10,000 flips

What is the right answer to the “coin question”?

theta <- (1:9)/10
prior <- rep(1/9, 9)
likelihood <- dbinom(6, size=10, prob=theta)
posterior <- prior*likelihood/sum(prior*likelihood)
names(posterior) <- theta
barplot(posterior, xlab='true prob of heads', main='posterior probability')

Motivating example

Now suppose we want to estimate the probability of heads for a coin without knowing the possible values. (or, a disease incidence, or error rate in an experiment, …)

We flip it $n$ times and get $z$ Heads.

The likelihood of this, given the prob-of-heads $\theta$, is: $$p(z|\theta )=(\genfrac{}{}{0ex}{}{n}{z}){\theta }^z(1-\theta {)}^{n-z}.$$

How to weight the possible $\theta$? Need a flexible set of weighting functions, i.e., prior distributions on $[0,1]$.

Beta-Binomial Bayesian analysis

If $$\begin{array}{rl}P& \sim \text{Beta}(a,b)\\ Z& \sim \text{Binom}(n,P),\end{array}$$

then “miraculously”,

$$\begin{array}{r}(P|Z=z)\sim \text{Beta}(a+z,b+n-z).\end{array}$$

Discuss:

We flip an odd-looking coin 100 times, and get 65 heads. What is it’s true* probability of heads?

1. What prior to use?

2. Plot the prior and the posterior.

3. Is it reasonable that $\theta =1/2$?

4. Best guess at $\theta$?

5. How far off are we, probably?

Tools include: rbeta( )

IN CLASS

prior_alpha <- prior_beta <- 1
post_alpha <- prior_alpha + 65
post_beta <- prior_beta + 35
post_mean <- post_alpha / (post_alpha + post_beta)

xvals <- seq(0, 1, length.out=101)
plot(xvals, dbeta(xvals, prior_alpha, prior_beta),
     type='l', ylim=c(0,10), ylab='density', xlab='theta')
lines(xvals, dbeta(xvals, post_alpha, post_beta),
    col='red')
abline(v=post_mean, lty=1, col='green')
abline(v= qbeta(c(.025, .975), post_alpha, post_beta),
       lty=3, col='blue')
legend("topleft", lty=c(1,1,1,3), col=c('black', 'red', 'green', 'blue'),
        legend=c("prior", "posterior", "posterior mean", "95% CI"))
abline(v=1/2, lty=3)

It doesn’t look likely that $\theta =1/2$: in fact, the posterior probability that $\theta \le 1/2$ is pbeta(1/2, post_alpha, post_beta). The posterior mean is 0.6470588. Our guess is probably acurate to within about 10% or so; a 95% credible interval is from 0.7363926 to 0.5522616.

If we were being frequentist, a 95% confidence interval for the mean (i.e., the proportion of heads) would be almost exactly the same thing.

phat <- 65/100
se <- sqrt(phat * (1 - phat) / 100)
conf_int <- phat + c(+1, -1) * 1.96 * se