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{aligned} \P\{ \text{6 H in 10 flips} \given \text{coin A} \} &= \binom{10}{6} (.75)^6 (.25)^4 \\ &= 0.146 \end{aligned}\] and \[\begin{aligned} \P\{ \text{6 H in 10 flips} \given \text{coin B} \} &= \binom{10}{6} (.25)^6 (.75)^4 \\ &= 0.016 \end{aligned}\] … so, probably coin (A)?

For a precise answer…

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

   \[\begin{aligned} \P\{ \text{coin A} \given \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{aligned}\]

Probability rules:

0. Probabilities are proportions: \(\hspace{2em} 0 \le \P\{A\} \le 1\)

1. Everything: \(\hspace{2em} \P\{ \Omega \} = 1\)

2. Complements: \(\hspace{2em} \P\{ \text{not } A\} = 1 - \P\{A\}\)

3. Disjoint events: If \(\hspace{2em} \P\{A \text{ and } B\} = 0\) then \(\hspace{2em} \P\{A \text{ or } B\} = \P\{A\} + \P\{B\}\).

4. Independence: \(A\) and \(B\) are independent iff \(\P\{A \text{ and } B\} = \P\{A\} \P\{B\}\).

5. Conditional probability: \[\P\{A \given B\} = \frac{\P\{A \text{ and } B\}}{ \P\{B\} }\]

Bayes’ rule

A consequence is

\[\P\{B \given A\} = \frac{\P\{B\} \P\{A \given B\}}{ \P\{A\} } .\]

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, \ldots, 0.8, 0.9,\) find \[\begin{aligned} \P\{\text{ coin has prob $\theta$ } \given \text{ 6 H in 10 flips } \} . \end{aligned}\]

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

The Beta Distribution

If \[P \sim \text{Beta}(a,b)\] then \(P\) has probability density \[p(\theta) = \frac{ \theta^{a-1} (1 - \theta)^{b-1} }{ B(a,b) } . \]

• Takes values between 0 and 1.

• If \(U_{(1)} < U_{(2)} < \cdots < U_{(n)}\) are sorted, independent \(\text{Unif}[0,1]\) then \(U_{(k)} \sim \text{Beta}(k, n-k+1)\).

• Mean: \(a/(a+b)\).

• Larger \(a+b\) is more tightly concentrated (like \(1/\sqrt{a+b}\))

Exercise

In pairs, simulate:

1. One thousand “random coins” whose probabilities are drawn from a \(\Beta(5,5)\) distribution. (rbeta()) Make a histogram of these probabilities.

2. Flip each coin ten times and record the number of heads. (rbinom())

3. Make a histogram of the probabilities of those coins that got exactly 3 heads, and compare to the first histogram.

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 \given \theta) = \binom{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]\).

What would we use if:

• the coin is probably close to fair.

• the disease is probably quite rare.

• no idea whatsoever.

Beta-Binomial Bayesian analysis

If \[\begin{aligned} P &\sim \text{Beta}(a,b) \\ Z &\sim \text{Binom}(n,P) , \end{aligned}\] then by Bayes’ rule: \[\begin{aligned} \P\{ P = \theta \given Z = z\} &= \frac{\P\{Z = z \given P = \theta \} \P\{P = \theta\}}{\P\{Z = z\}} \\ &= \frac{ \binom{n}{z}\theta^z (1-\theta)^{n-z} \times \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)} }{ \text{(something)} } \\ &= \text{(something else)} \times \theta^{a + z - 1} (1-\theta)^{b + n - z - 1} . \end{aligned}\]

Discuss/demonstration:

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

1. True = ??

2. What prior to use?

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

4. Best guess at \(\theta\)?

5. How far off are we, probably?

How do we communicate results?

If we want a point estimate:

1. posterior mean,
2. posterior median, or
3. maximum a posteriori estimate (“MAP”: highest posterior density).

These all convey “where the posterior distribution is”, more or less.

Credible region

Definition: A 95% credible region is a portion of parameter space having a total of 95% of the posterior probability.

Interpretation #1

If we construct a 95% credible interval for \(\theta\) for each of many datasets; and the coin in each dataset has \(\theta\) drawn independently from the prior, then the true \(\theta\) will fall in its credible interval for 95% of the datasets.

Interpretation #2

If we construct a 95% credible interval for \(\theta\) with a dataset, and the distribution of the coin’s true \(\theta\) across many parallel universes is given by the prior, then the true \(\theta\) will be in the credible interval in 95% of those universes.

Interpretation #3

Given my prior beliefs (prior distribution), the posterior distribution is the most rational\({}^*\) way to update my beliefs to account for the data.

But which credible interval?

Definition: The “95% highest density interval” is the 95% credible interval with the highest posterior probability density at each point.

Motivating problem: more coins

Suppose now we have data from \(n\) different coins from the same source. We don’t assume they have the same \(\theta\), but don’t know what it’s distribution is, so try to learn it.

\[\begin{aligned} Z_i &\sim \Binom(N_i, \theta_i) \\ \theta_i &\sim \Beta(\text{mode}=\omega, \text{conc}=\kappa) \\ \omega & = ? \\ \kappa & = ? \end{aligned}\]

note: The “mode” and “concentration” are related to the shape parameters by: \(\alpha = \omega (\kappa - 2) + 1\) and \(\beta = (1 - \omega) (\kappa - 2) + 1\).

Binomial versus Beta-Binomial

What is different between:

1. Pick a value of \(\theta\) at random from \(\Beta(3,1)\). Flip one thousand \(\theta\)-coins, 500 times each.

2. Pick one thousand random \(\theta_i \sim \Beta(3,1)\) values. Flip one thousand coins, one for each \(\theta_i\), 500 times each.

Stan

The skeletal Stan program

data {
    // stuff you input
}
transformed data {
    // stuff that's calculated from the data (just once, at the start)
}
parameters {
    // stuff you want to learn about
}
transformed parameters {
    // stuff that's calculated from the parameters (at every step)
}
model {
    // the action!
}
generated quantities {
    // stuff you want computed also along the way
}

First, in words:

We’ve flipped a coin 10 times and got 6 Heads. We think the coin is close to fair, so put a \(\Beta(20,20)\) prior on it’s probability of heads, and want the posterior distribution.

\[\begin{aligned} Z &\sim \Binom(10, \theta) \\ \theta &\sim \Beta(20, 20) \end{aligned}\] What’s our best guess at \(\theta\)?

The Stan model

library(rstan)
stan_block <- "
data {
    int N;   // number of flips
    int Z;   // number of heads
}
parameters {
    // probability of heads
    real<lower=0,upper=1> theta;
}
model {
    Z ~ binomial(N, theta);
    theta ~ beta(20, 20);
}
"
bb_model <- stan_model(
               model_code=stan_block)

Optimization: maximum likelihood

With a uniform prior, the “maximum posterior” parameter values are also the maximum likelihood values.

> help(optimizing)

optimizing                package:rstan                R Documentation

Obtain a point estimate by maximizing the joint posterior

Description:

     Obtain a point estimate by maximizing the joint posterior from the
     model defined by class 'stanmodel'.

Usage:

     ## S4 method for signature 'stanmodel'
     optimizing(object, data = list(),

Maximum posterior

library(rstan)
(fit <- optimizing(bb_model,  # stan model from above
                   data=list(N=10, Z=6)))

## Initial log joint probability = -42.3369
## Optimization terminated normally: 
##   Convergence detected: relative gradient magnitude is below tolerance

## $par
##     theta 
## 0.5208333 
## 
## $value
## [1] -33.22939
## 
## $return_code
## [1] 0

The answer!

post_fun <- function (p) {
    n <- 10; z <- 6; a <- b <- 20
    lik <- dbinom(z, size=n, prob=p)
    prior <- dbeta(p, a, b)
    return( prior * lik )
}
curve(post_fun, 0, 1, xlab=expression(theta), ylab='posterior prob')
points(fit$par, post_fun(fit$par), col='red', pch=20)

Summary

posterior \(=\) prior \(\times\) likelihood

Exercise: check the Beta-Binomial

1. Simulate \(10^6\) coin probabilities, called \(\theta\), from Beta(5,5). (rbeta)

2. For each coin, simulate 10 flips. (rbinom)

3. Make a histogram of the true probabilities (values of \(\theta\)) of only those coins having 3 of 10 heads.

4. Compare the distribution to Beta(\(a\),\(b\)) – with what \(a\) and \(b\)? (dbeta)

5. Explain what happened.