Bayesian inference by conditional simulation

Today: Bayesian statistics intro

So far, we have

• Formulated generative models for observed data
• Found most plausible parameter values using SMM
• Looked at the uncertainty associated with the parameter estimates using the parametric bootstrap

Bayesian statistics provides a different way of assessing plausible parameter values given a generative model.

The Bayesian idea

We are still going to have a data model: \[ X \sim P_\theta \] for some unknown \(\theta \in \Theta\).

The new thing is that we augment this process: we imagine that \(\theta\) is also a random variable that was drawn from some distribution.

\[ \theta \sim \pi \]

We call \(\pi\) the “prior.”

How does this help?

We have specified the joint distribution of \((\theta, X)\).

Therefore, we can obtain the conditional distribution of \(\theta \mid X = x_{\text{obs}}\).

This is called the “posterior” distribution. It describes how plausible various values of \(\theta\) are after observing \(X=x_{\text{obs}}\), under the model we specified.

We already know how to sample from conditional distributions

In the disease-test example, we simulated disease status \(D\) and test result \(T\), then kept runs with \((D,T)\in \{(0,1), (1,1)\}\), the outcomes with a positive test, and examined their \(D\) values.

For Bayesian inference:

• \(\theta\) takes the place of the unknown quantity \(D\)
• the observed data take the place of \(T\).

When \(X=x_{\text{obs}}\) has positive probability, we can sample from the posterior distribution by rejection:

• Draw \(\theta \sim \pi\), draw \(X \sim P_\theta\)
• Keep runs with \(X = x_{\text{obs}}\)

The values of \(\theta\) that we keep will be draws from the posterior distribution.

The distribution of many retained values approximates the posterior.

Example 1: Three possible coins

Suppose we are able to observe two flips of a coin.

Suppose further that we know that the coin’s probability of heads \(\theta\) is one of \(\{\frac14, \frac12, \frac34 \}\), and that each of those options is equally likely (maybe we have this information because we are playing some carnival game).

We then flip the coin twice and observe HH.

After we’ve seen HH, how plausible are the various possible values for \(\theta\)?

Write down the generative model

First the prior: \[ \theta\sim\pi, \qquad \pi(\{1/4\}) = \pi(\{1/2\}) = \pi(\{3/4\}) = 1/3 \]

Then the probability of the data given \(\theta\): \[ X_1, X_2\stackrel{\mathrm{iid}}{\sim} \text{Bernoulli}(\theta), \] with \(X_i = 0\) denoting tails on flip \(i\) and \(X_i = 1\) denoting heads on flip \(i\).

Rejection sampling for posterior distributions

One run of the enlarged generative model is:

1. Draw \(\theta \sim \pi\) (uniformly from \(\{1/4,1/2,3/4\}\)).
2. Using that value of \(\theta\), flip the coin twice independently.
3. Record the outcome \((\theta,X_1,X_2)\).
4. If \(X_1 = 1, X_2 = 1\), then record \(\theta\). Otherwise return to 1 and try again until we get an acceptance.

Three-coin example, step 1 in code

set.seed(681)
simulate_three_coin <- function(prior_probs = c(1/3, 1/3, 1/3)) {
  theta <- sample(c(0.25, 0.5, 0.75),
                  size = 1, prob = prior_probs)
  X <- rbinom(2, size = 1,
              prob = theta)
  return(list(theta = theta, X = X))
}

import numpy as np
rng = np.random.default_rng(681)

def simulate_three_coin(prior_probs=(1/3, 1/3, 1/3)):
    values = [0.25, 0.5, 0.75]
    theta = rng.choice(values, p=prior_probs)
    X = rng.binomial(1, theta, 2)
    return {"theta": theta, "X": X}

Three-coin example, steps 2-4: Simulate from unconditional, check condition, repeat until an acceptance

rejection_sample_three_coin <-
  function(prior_probs = c(1/3, 1/3, 1/3)) {
    while (TRUE) {
      run <- simulate_three_coin(prior_probs)
      if (all(run$X == c(1, 1)))
        return(run$theta)
    }
  }
rejection_sample_three_coin()

def rejection_sample_three_coin(prior_probs=(1/3, 1/3, 1/3)):
    while True:
        run = simulate_three_coin(prior_probs)
        if list(run["X"]) == [1, 1]:
            return run["theta"]

rejection_sample_three_coin()

Three-coin example: estimating the posterior

draws <- replicate(5000,
  rejection_sample_three_coin())

# Fraction of accepted draws
# at each value:
table(draws) / 5000

draws = np.array(
    [rejection_sample_three_coin()
     for _ in range(5000)])

# Fraction of accepted draws
# at each value:
for theta in [0.25, 0.5, 0.75]:
    frac = np.mean(draws == theta)
    print(theta, frac)

The fraction at each value estimates its posterior probability. With more accepted draws, these estimates become more stable.

What if the prior had been different?

Suppose our prior probabilities for \(\theta\in\{1/4,1/2,3/4\}\) were instead \((0.6,0.3,0.1)\).

Use those probabilities when drawing \(\theta\), then run the same rejection sampler. The coin-flip model and the observed HH stay the same.

new_prior <- c(0.6, 0.3, 0.1)
new_draws <- replicate(5000,
  rejection_sample_three_coin(new_prior))
table(factor(new_draws,
  levels = c(0.25, 0.5, 0.75))) / 5000

new_prior = [0.6, 0.3, 0.1]
new_draws = np.array([
    rejection_sample_three_coin(new_prior)
    for _ in range(5000)])
for theta in [0.25, 0.5, 0.75]:
    print(theta, np.mean(new_draws == theta))

Compare these frequencies with the equal-prior results. Now \(\theta=1/2\) is typically the most frequent value among the accepted draws.

What have we gained?

We can now turn an unknown parameter into part of a generative model:

\[ \text{draw parameter}\quad\longrightarrow\quad \text{simulate data given parameter}. \]

Conditioning on \(X = x_{\text{obs}}\), the outcomes whose data match the observation, gives a posterior distribution over the parameter.

We estimated that distribution by retaining simulated runs. But requiring \(\theta\) to be one of only three values was artificial. What if any value between 0 and 1 is possible?

Example 2: A prior with infinitely many values

Keep the same coin-flipping story, but now draw

\[ \theta\sim\operatorname{Uniform}(0,1). \]

The interpretation here is that this prior is supposed to convey that before we see the data, we have no particular belief about what \(\theta\) might be.

Suppose we flip the coin six times and observe the complete sequence HHTHHT (four heads and two tails).

The parameter now has infinitely many possible values, but we can still simulate \(\theta\) and the data together.

Rejection sampling still works

Let \(X=(X_1,\dots,X_6)\) record the six flips, with heads coded as 1. The observed vector is \(x_{\mathrm{obs}}=(1,1,0,1,1,0)\). One run of the joint model gives:

\[ \theta\sim\operatorname{Uniform}(0,1), \qquad X_1,\dots,X_6\stackrel{\mathrm{iid}}{\sim}\operatorname{Bernoulli}(\theta), \]

where the flips use the value of \(\theta\) drawn first. Simulate many independent runs, retain only those with \(X=x_{\mathrm{obs}}\), and inspect their \(\theta\) values.

Because \(X\) is discrete, \(\Pr(X=x_{\mathrm{obs}})>0\). This is exact conditioning on all six observed flips.

Continuous-prior coin, step 1 in code

set.seed(681)
simulate_continuous_coin <- function() {
  theta <- runif(1)
  X <- rbinom(6, size = 1,
              prob = theta)
  return(list(theta = theta, X = X))
}

import numpy as np
rng = np.random.default_rng(681)

def simulate_continuous_coin():
    theta = rng.uniform()
    X = rng.binomial(1, theta, 6)
    return {"theta": theta, "X": X}

Continuous-prior coin, steps 2-4: Simulate, check, repeat until acceptance

rejection_sample_continuous_coin <-
  function() {
    while (TRUE) {
      run <- simulate_continuous_coin()
      if (all(run$X == c(1, 1, 0, 1, 1, 0)))
        return(run$theta)
    }
  }
rejection_sample_continuous_coin()

def rejection_sample_continuous_coin():
    while True:
        run = simulate_continuous_coin()
        if list(run["X"]) == [1, 1, 0, 1, 1, 0]:
            return run["theta"]

rejection_sample_continuous_coin()

Each accepted \(\theta\) is a draw from the posterior distribution of \(\theta\) given \(X=x_{\mathrm{obs}}\).

Continuous-prior coin: estimating the posterior

draws <- replicate(2000,
  rejection_sample_continuous_coin())
mean(draws)
hist(draws, breaks = 40,
     xlim = c(0, 1),
     xlab = "theta", main = "")

import matplotlib.pyplot as plt
draws = np.array(
    [rejection_sample_continuous_coin()
     for _ in range(2000)])
print(draws.mean())
plt.hist(draws, bins=40)
plt.xlim(0, 1)
plt.xlabel("theta")
plt.show()

Compare prior and posterior

The retained values favor larger \(\theta\), while values across the interval remain possible.

Interpretation of the retained samples

• The samples we obtain through rejection sampling are supposed to be draws from the posterior
• We can use the fraction of samples that satisfy various properties as estimates of the probability a draw from the posterior satisfies those properties (e.g. \(P(\theta \in [.1,.2] \mid X = x_{\text{obs}})\) is approximated by the fraction of retained \(\theta\) values fall in \([.1, .2]\)).
• We usually interpret the distribution as giving plausibilities to various values of \(\theta\).

Can we be more efficient?

Even at the most plausible value of \(\theta\), any particular sequence of \(X_i\) values is going to have a small probability, and rejection sampling can be inefficient.

Can we get away with conditioning on a set with a larger probability?

Let \(S=\sum_{i=1}^6 X_i\) be the number of heads. Our observation has \(S=4\). We now retain every simulated sequence with four heads.

rejection_sample_coin_count <-
  function() {
    while (TRUE) {
      run <- simulate_continuous_coin()
      if (sum(run$X) == 4)
        return(run$theta)
    }
  }
draws_count <- replicate(2000,
  rejection_sample_coin_count())

def rejection_sample_coin_count():
    while True:
        run = simulate_continuous_coin()
        if sum(run["X"]) == 4:
            return run["theta"]

draws_count = np.array(
    [rejection_sample_coin_count()
     for _ in range(2000)])

The two conditions give the same posterior

The same comparison with a Q-Q plot

Sufficient statistics

• For independent flips with a common head probability, conditioning on \(S=4\) gives the same distribution of \(\theta\) as conditioning on the complete sequence HHTHHT.
• Once we know there were four heads, every possible order of those heads is equally likely, whatever \(\theta\) is. The particular order adds no information about \(\theta\).
• We call \(S\) a sufficient statistic for \(\theta\) in this model.
• Formal definition (optional): A statistic \(S(X)\) is sufficient for \(\theta\) if the conditional distribution of the full data \(X\), given \(S(X)\), does not depend on \(\theta\).

The simpler condition accepts more runs

There are \({6\choose4}=15\) sequences with four heads. Each has the same probability for any fixed \(\theta\).

Matching the count therefore accepts 15 times as many runs as matching HHTHHT, while targeting the same posterior distribution.

Example 3: Rainy and dry days

Let \(X_d=1\) mean rain on day \(d\) and \(X_d=0\) mean dry weather. For six consecutive days, use the model

\[ \begin{aligned} \theta &\sim \operatorname{Uniform}(0,1), \\ X_1\mid\theta &\sim \operatorname{Bernoulli}(\theta), \\ Z_d &\overset{\mathrm{iid}}{\sim} \operatorname{Bernoulli}(1/2), \qquad d=2,\ldots,6, \\ X_d\mid X_{d-1},Z_d,\theta &\sim \operatorname{Bernoulli}\!\left( Z_dX_{d-1}+(1-Z_d)\theta\right), \qquad d=2,\ldots,6. \end{aligned} \]

Each \(Z_d\) is drawn independently of \(\theta\) and earlier days. If \(Z_d=1\), then \(X_d=X_{d-1}\); if \(Z_d=0\), then \(X_d\) is a fresh Bernoulli\((\theta)\) draw.

The persistence probability is known. Only \(\theta\) is unknown, and each day has marginal rain probability \(\theta\) even though consecutive days are dependent.

The order of rainy days matters

Here \(R\) means rainy and \(D\) means dry. Both DRRRRD and RDRDRR have four rainy days out of six.

In DRRRRD, one rainy spell begins on day 2. The next three rainy days could result from persistence.

In RDRDRR, rain begins on days 1, 3, and 5. Each return after a dry day requires a fresh rainy draw.

The three-start sequence therefore gives more evidence for a high \(\theta\) than the one-start sequence, despite the identical rain count.

Rain example, step 1 in code: R

As for the coin, give \(\theta\) a Uniform\((0,1)\) prior. One run draws \(\theta\), then six days using that \(\theta\). Write \(X=(X_1,\dots,X_6)\), with \(X_d=1\) if day \(d\) is rainy.

set.seed(681)
simulate_rain <- function() {
  theta <- runif(1)
  X <- integer(6)
  X[1] <- rbinom(1, size = 1, prob = theta)
  for (day in 2:6) {
    if (runif(1) < 0.5) {
      X[day] <- X[day - 1]                         # weather persists
    } else {
      X[day] <- rbinom(1, size = 1, prob = theta)   # new state
    }
  }
  return(list(theta = theta, X = X))
}

Rain example, step 1 in code: Python

import numpy as np
rng = np.random.default_rng(681)

def simulate_rain():
    theta = rng.uniform()
    X = [rng.binomial(1, theta)]
    for day in range(1, 6):
        if rng.uniform() < 0.5:
            X.append(X[-1])                    # weather persists
        else:
            X.append(rng.binomial(1, theta))  # new state
    return {"theta": theta, "X": X}

Rain example, steps 2-4: Simulate from unconditional, check condition, repeat until an acceptance

rejection_sample_rain <-
  function(E) {
    while (TRUE) {
      run <- simulate_rain()
      if (E(run$X))
        return(run$theta)
    }
  }

def rejection_sample_rain(E):
    while True:
        run = simulate_rain()
        if E(run["X"]):
            return run["theta"]

Here E is a function that checks whether the simulated days are in the set we condition on. The same sampler works for any choice of \(E\).

Rain example: matching each complete sequence

E_DRRRRD <- function(X)
  all(X == c(0, 1, 1, 1, 1, 0))
E_RDRDRR <- function(X)
  all(X == c(1, 0, 1, 0, 1, 1))

rejection_sample_rain(E_RDRDRR)

def E_DRRRRD(X):
    return X == [0, 1, 1, 1, 1, 0]

def E_RDRDRR(X):
    return X == [1, 0, 1, 0, 1, 1]

rejection_sample_rain(E_RDRDRR)

\(E_{\text{DRRRRD}}\) and \(E_{\text{RDRDRR}}\) each match one complete six-day sequence.

Rain example: comparing the exact posteriors

draws <- function(E)
  replicate(2000,
    rejection_sample_rain(E))

mean(draws(E_DRRRRD))
mean(draws(E_RDRDRR))

def draws(E):
    return np.array([
        rejection_sample_rain(E)
        for _ in range(2000)])

print(draws(E_DRRRRD).mean())
print(draws(E_RDRDRR).mean())

The simulated posterior mean is higher after RDRDRR than after DRRRRD. The estimates vary from run to run because they use finitely many accepted draws.

Posteriors for two weather patterns

Both sequences contain four rainy days, but their posteriors differ. The RDRDRR posterior is shifted toward larger \(\theta\) values than the DRRRRD posterior.

Can we get the posterior by conditioning on the number of rainy days?

• We saw above that in both cases, \(S = \sum_{d=1}^6 X_d = 4\), but they gave different posteriors.
• If we condition only on \(S=4\), we cannot recover the posterior given either complete sequence. (Technically, because \(S\) is not sufficient for \(\theta\) in this model.)

What happens if we condition on the count anyway?

E_count <- function(X)
  sum(X) == 4
count_draws <- draws(E_count)
mean(count_draws)

def E_count(X):
    return sum(X) == 4

count_draws = draws(E_count)
print(count_draws.mean())

Posterior given four rainy days, ignoring order

Summing up

• Rejection sampling where you match exactly the observed data gets very costly
• In simple models, you can condition on a sufficient statistic instead and get more acceptances, but not all models are simple enough for you to be able to compute one

What we can now do

For a generative model with an unknown parameter:

1. Specify a prior distribution for the parameter.
2. Simulate the parameter and a dataset from the joint model.
3. Retain runs with \(X = x_{\text{obs}}\).
4. Use the retained parameters to study the posterior distribution.

The three-coin example let us estimate posterior probabilities by counting retained values. With a continuous prior, the same conditional-simulation procedure gives a distribution of retained parameter values.

What is still a problem?

• Rejection sampling relies on the probability of the observed data being positive.
• If we have a continuous observation model, probability of the observed data will be zero.
• Even with discrete data, if we have a lot of it, acceptance probabilities will be low and rejection sampling will become infeasible.
• Sufficient statistics can help us in some cases, but not always
• ABC (next time) will help with some of this