Conditional distributions and simulation

Today

Conditional distributions and simulation from conditional distributions via rejection sampling

Recall: Random variables and random vectors

A random variable represents a numerical outcome that is not yet known. We use a capital letter, such as \(X\), for the random variable and a lowercase letter, such as \(x\), for a possible value.

A random vector collects several random variables, such as \((X,Y)\).

For example, \(X\) might indicate whether a person has a disease, and \(Y\) might indicate whether their test is positive.

The pair \((X,Y)\) records both outcomes for the same person.

Recall: Probability laws/distributions

A probability law, or distribution, specifies the probabilities of events involving a random variable or random vector.

For a single random variable, these might be probabilities such as

\[ P(X=1) \qquad \text{or} \qquad P(a<X\le b). \]

For a random vector, we call its distribution a joint distribution. It specifies probabilities of simultaneous statements about its components, such as

\[ P(X=1,\ Y=1). \]

Recall: Marginal distributions

The distribution of one component of a random vector considered by itself is its marginal distribution.

The joint distribution of \((X,Y)\) tells us how likely different combinations of \(X\) and \(Y\) are.

The marginal distribution of \(X\) tells us how likely different values of \(X\) are, regardless of the value of \(Y\).

To simulate from this marginal distribution, generate \((X,Y)\) jointly and keep only \(X\) from every run.

Conditional distributions: why?

Suppose we have partial information about the outcome of a random variable.

How can we use that to update our understanding of what values the random variable can take?

“Partial information” is interpreted as “we know the random variable/random vector’s value falls in a certain set.”

Conditional distributions: definition

Suppose \(X\) is a discrete random variable (or or discrete random vector), \(E\) is some subset of the values \(X\) can take, and \(\Pr(X \in E) > 0\).

Then the conditional distribution of \(X\) given \(X \in E\) is a new probability distribution, denoted \(X \mid X \in E\) (the distribution of \(X\) given \(X \in E\)).

The probabilities it assigns to possible outcomes are

\[\begin{align*} \Pr(X=x\mid X \in E) &=\frac{\Pr(X=x\text{ and }X \in E)}{\Pr(X \in E)} \\ &= \begin{cases} \frac{\Pr(X = x)}{ \Pr(X \in E)} & x \in E \\ 0 & x \not \in E \end{cases}. \end{align*}\]

Notes

• The idea is that we are making a new probability distribution where the only possible values are the elements of \(E\).
• The probabilities given by the conditional distribution are then the original probabilities for \(X\) on the restricted collection renormalized so that they sum to 1.
• It is probably more common to say that \(E\) is some “event,” e.g. \(E\) is the event that \(X\) falls in the set \(S\), and then the conditional distribution is written \(X \mid E\), and we say \(\Pr(X = x \mid E) = \Pr(X = x \text{ and }E) / \Pr(E)\).
• For continuous random variables, the ideas are the same but the math is more fiddly.

Example 1: How much longer will symptoms last?

Suppose we are modeling the cessation of the nausea and vomiting of pregnancy.

We will formalize this by saying that the random variable \(R\) describes the week of cessation of symptoms, with allowable values and probabilities given in the table below

Numbers from: Gadsby, Barnie-Adshead & Jagger (1993), Table 1.

Suppose that the woman is experiencing symptoms at the beginning of week 14.

What is the chance her symptoms will end before the beginning of week 15?

Example 2: Probability of having a disease given positive test

Suppose we have a test for a disease. In this hypothetical example:

• 1% of people in the population have the disease.
• Among people with the disease, 90% test positive.
• Among people without the disease, 5% test positive.

We select a person at random from this population and they test positive.

What is the probability that they have the disease?

Let \(D\) indicate whether the randomly selected person has the disease:

• \(D=1\) if they have the disease.
• \(D=0\) if they do not.

Let \(T\) be the result of the test:

• \(T = 1\) if the test is positive,
• \(T = 0\) if the test is negative.

The previous page allows us to write down the joint distribution of the random vector \((D,T)\). The probabilities of all four combinations are:

Example 3: A series that is already underway

Two teams play a best-of-five series. The first to win three games wins the series.

The first three games have been played:

• Team A won game 1.
• Team B won game 2.
• Team A won game 3.

What is the probability that Team A wins the series, given these results?

Let \(X_i\) indicate the outcome of game \(i\):

• \(X_i=1\) if Team A wins.
• \(X_i=0\) if Team B wins.

Let \(W\) indicate the series winner:

• \(W=1\) if Team A wins the series.
• \(W=0\) if Team B wins the series.

We describe the complete outcome using the random vector

\[ \mathbf{X}=(X_1,X_2,X_3,X_4,X_5,W). \]

For convenience, we include outcomes for all five games, even if the series ends earlier. Outcomes after one team has won three games do not affect the series winner.

Our model assigns probability \(1/32\) to each of the 32 possible five-game sequences.

Each sequence determines the winner: \(W=1\) if at least three game outcomes are \(1\), and \(W=0\) otherwise.

Thus, the augmented vector \(\mathbf{X}\) has 32 possible values, each with probability \(1/32\). For example:

The last row is impossible: the game outcomes say A wins, but the recorded series winner says B wins.

Let \(E\) be the set of possible complete outcomes consistent with the three games we observed:

\[ E=\{ (1,0,1,0,0,0),\ (1,0,1,0,1,1),\ (1,0,1,1,0,1),\ (1,0,1,1,1,1) \}. \]

Knowing the results of the first three games tells us that \(\mathbf{X}\in E\).

We want the conditional distribution of the winner,

\[ W\mid\mathbf{X}\in E. \]

Each of the four vectors in \(E\) has probability \(1/32\) under the joint distribution.

Adding these probabilities gives

\[ P(\mathbf{X}\in E)=\frac{4}{32}. \]

There are three rows with \(W=1\) and one row with \(W=0\).

To obtain the conditional distribution, we restrict to outcomes in \(E\) and divide by their total probability:

\[ P(W=1\mid\mathbf{X}\in E) =\frac{3/32}{4/32} =\frac34, \]

\[ P(W=0\mid\mathbf{X}\in E) =\frac{1/32}{4/32} =\frac14. \]

Given the observed results A, B, A, Team A has a 75% chance of winning the series under this model.

Rejection sampling to simulate from a conditional distribution

Suppose we can simulate complete outcomes of \(X\), where \(X\) may be a random vector.

To draw from the conditional distribution of \(X\) given \(X \in E\), where \(\Pr(X \in E) > 0\):

1. Run the simulator.
2. If \(X \in E\), keep the value of \(X\) from that run.
3. Otherwise, discard the run.
4. Repeat independently until we get an acceptance.

This is rejection sampling. The accepted value will come from the desired conditional distribution.

Correspondence between rejection sampling and definition

To obtain a conditional distribution, we restrict to outcomes in \(E\) and rescale their probabilities to sum to one.

Rejection sampling implements this by simulating from the full distribution and keeping only simulated outcomes in \(E\).

What do we retain in each example?

We always examine the value of interest from the same run whose outcome was in \(E\).

Let’s implement this for the disease-testing example.

Disease example, step 1: Simulate from unconditional distribution

We could use the joint distribution described above directly, or we could write the joint as the following model:

\[\begin{align*} D&\sim\operatorname{Bernoulli}(0.01) \\ T&\sim \begin{cases} \operatorname{Bernoulli}(0.90), & \text{if } D=1,\\ \operatorname{Bernoulli}(0.05), & \text{if } D=0. \end{cases} \end{align*}\]

Disease example, step 1 in code

set.seed(681)
simulate_DT <- function() {
  D <- rbinom(1, size = 1,
              prob = 0.01)
  p <- ifelse(D == 1,
              0.90, 0.05)
  T <- rbinom(1, size = 1,
              prob = p)
  return(list(D = D, T = T))
}

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

def simulate_DT():
    D = rng.binomial(1, 0.01)
    p = 0.90 if D == 1 else 0.05
    T = rng.binomial(1, p)
    return {"D": D, "T": T}

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

rejection_sample_DT <-
  function() {
    while (TRUE) {
      DT <- simulate_DT()
      if (DT$T == 1)
        return(DT$D)
    }
  }
rejection_sample_DT()

def rejection_sample_DT():
    while True:
        DT = simulate_DT()
        if DT["T"] == 1:
            return DT["D"]

rejection_sample_DT()

Duration example, step 1: Simulate from unconditional distribution

For the duration example, one run produces the cessation week \(R\) alone. Its distribution is the table from Example 1:

We keep runs with \(R \in E\), where \(E = \{14, 15, 16\}\).

Duration example, step 1 in code

set.seed(681)
simulate_R <- function() {
  probs <- c(0.082, 0.092, 0.099,
             0.151, 0.175, 0.086,
             0.110, 0.072, 0.041,
             0.092)
  R <- sample(7:16, size = 1,
              prob = probs)
  return(R)
}

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

def simulate_R():
    probs = [0.082, 0.092, 0.099,
             0.151, 0.175, 0.086,
             0.110, 0.072, 0.041,
             0.092]
    weeks = np.arange(7, 17)
    R = rng.choice(weeks, p=probs)
    return R

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

rejection_sample_R <-
  function() {
    while (TRUE) {
      R <- simulate_R()
      if (R %in% c(14, 15, 16))
        return(R)
    }
  }
rejection_sample_R()

def rejection_sample_R():
    while True:
        R = simulate_R()
        if R in [14, 15, 16]:
            return R

rejection_sample_R()

Series example, step 1: Simulate from unconditional distribution

For the series example, one run plays all five games and records the winner, giving \(\mathbf{X}=(X_1,\dots,X_5,W)\):

\[\begin{align*} X_1, \dots, X_5 &\stackrel{\text{iid}}{\sim} \operatorname{Bernoulli}(1/2) \\ W &= \begin{cases} 1, & \text{if } X_1 + \dots + X_5 \geq 3,\\ 0, & \text{otherwise.} \end{cases} \end{align*}\]

We keep runs with \(\mathbf{X} \in E\): the first three games are A, B, A, that is, \((X_1, X_2, X_3) = (1, 0, 1)\).

Series example, step 1 in code

set.seed(681)
simulate_series <- function() {
  X <- rbinom(5, size = 1,
              prob = 0.5)
  W <- as.integer(sum(X) >= 3)
  return(list(X = X, W = W))
}

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

def simulate_series():
    X = rng.binomial(1, 0.5, 5)
    W = int(X.sum() >= 3)
    return {"X": X, "W": W}

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

rejection_sample_series <-
  function() {
    while (TRUE) {
      XW <- simulate_series()
      start <- XW$X[1:3]
      if (all(start == c(1,0,1)))
        return(XW$W)
    }
  }
rejection_sample_series()

def rejection_sample_series():
    while True:
        XW = simulate_series()
        start = list(XW["X"][:3])
        if start == [1, 0, 1]:
            return XW["W"]

rejection_sample_series()

What we saw

• Conditioning on \(X \in E\) restricts attention to the outcomes in \(E\) and renormalizes probabilities within it.
• Rejection sampling implements that restriction by retaining runs whose outcome is in \(E\). We use the retained values to estimate conditional probabilities.
• The same procedure lets us interpret a test result, predict when symptoms will end, or update the chances of winning a series.

Next we will use the same idea for Bayesian inference: generate a parameter and data together, then condition on the outcome falling in the set \(E\) of outcomes whose data match what we observed.

Optional: why rejection sampling works

For students who have taken probability separately. We have not developed all the background needed to justify these steps in this class. You are not expected to reproduce this proof.

Run independent copies of the simulator. Write \(X_i\) for the outcome of run \(i\), let \(p=\Pr(X\in E)>0\), and let \(\tau\) be the first index with \(X_\tau\in E\). Since \(\Pr(\tau>n)=(1-p)^n\to0\), acceptance occurs almost surely.

For any set \(B\) of possible values of \(X\),

\[ \begin{aligned} \Pr(X_\tau\in B) &=\sum_{n=1}^{\infty}\Pr(\tau=n,\;X_n\in B)\\ &=\sum_{n=1}^{\infty}(1-p)^{n-1}\Pr(X\in B,\;X\in E)\\ &=\frac{\Pr(X\in B,\;X\in E)}{p} =\Pr(X\in B\mid X\in E). \end{aligned} \]

The first accepted draw therefore has exactly the desired conditional law.