Lab 4: Bayesian Inference

Course: INS-605: Data Analysis II
Lecturer: Sothea HAS, PhD

---

• Student name: …
• ID: …

---

• The jupyter notebook for this lab can be downloaded from here.

0. Before You Start

In lecture, we saw that Bayesian inference combines a prior belief \(p(\theta)\) with the likelihood of observed data \(p(x\mid\theta)\) to get a posterior \(p(\theta\mid x)\propto p(x\mid\theta)p(\theta)\).

In this lab, you will not need to derive any calculus by hand. Instead, we will:

• Recognize a real-world problem as a Beta-Binomial model (and, as a bonus, a Poisson-Gamma model).
• Use the conjugate update rules given below to jump straight from prior + data to posterior.
• Use scipy.stats and simulation (numpy.random) to visualize, summarize, and make decisions from the posterior.

TipBeta-Binomial Model

Conjugate update rule you will use everywhere in Section 1-3 (memorize this box!)

If \(\theta \sim \text{Beta}(\alpha, \beta)\) (prior) and we observe \(x\) successes out of \(n\) trials (data, Binomial likelihood), then:

\[(\theta \mid x) \sim \text{Beta}(\alpha + x,\ \beta + n - x) \quad \text{(posterior)}\]

That’s it — no integrals needed. \(\alpha\) and \(\beta\) act like “pseudo-counts” of successes/failures you already believed in before seeing the data.

Run the cell below to load the packages we need.

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import beta, gamma, binom, poisson

np.random.seed(605)
plt.rcParams['figure.figsize'] = (7, 4)
plt.rcParams['axes.grid'] = True

---

1. Problem Setting: Did the New Checkout Page Work?

An e-commerce company redesigned its checkout page and wants to know whether it improved the conversion rate \(\theta\) (the true, unknown proportion of visitors who complete a purchase).

• The old checkout page had a well-known conversion rate of about 20%.
• The company shows the new page to \(n = 50\) visitors, and \(x = 14\) of them purchase something.

We want to estimate \(\theta\), the true long-run conversion rate of the new page, and decide whether it is genuinely better than 20%.

1.1. Warm-up: The Frequentist (MLE) Estimate

A. Using the data above (\(n=50\), \(x=14\)), compute the sample proportion \(\hat\theta = x/n\). This is the Maximum Likelihood Estimate (MLE) of \(\theta\) (I dare you to prove this!).

B. In a markdown cell, briefly answer: with only 50 visitors, would you fully trust \(\hat\theta\) as “the” conversion rate? Why might it be risky to make a business decision based on this single number alone?

# A. Sample proportion (MLE)
n = 50
x = 14

theta_hat = x/n # TODO: replace None with it's actual value.

print(f"MLE estimate of theta: {theta_hat}")

MLE estimate of theta: 0.28

B. Your answer here:

…

---

2. Building the Bayesian Model: Prior, Likelihood, Posterior

Instead of relying only on 50 new visitors, let’s bring in what we already know: the old page converted at ~20%. We encode this belief as a prior on \(\theta\).

• Model: each visitor purchases independently with probability \(\theta\), so \(x \mid \theta \sim \text{Binomial}(n, \theta)\).
• Prior: we choose \(\theta \sim \text{Beta}(\alpha_0, \beta_0)\) with \(\alpha_0 = 4\), \(\beta_0 = 16\). What is the mean of this distribution? (Hint: the mean of a \(\text{Beta}(\alpha_0, \beta_0)\) is \(\frac{\alpha_0}{\alpha_0 + \beta_0}\).).

2.1. Visualize the Prior

A. Complete the code below to plot the density of \(\text{Beta}(4, 16)\) over \(\theta \in [0,1]\) using beta.pdf.

# A. Plot the prior Beta(4, 16)
alpha0, beta0 = 4, 16

# To do
theta0_hat = alpha0/(alpha0+beta0)

print(f'Mean of our prior belief: {theta0_hat}').

Mean of our prior belief: 0.2

• Is this a weak or strong prior? Why?

2.2. Update to the Posterior

B. Using the conjugate update rule from Section 0 (\(\alpha_{\text{post}} = \alpha_0 + x\), \(\beta_{\text{post}} = \beta_0 + n - x\)), compute the posterior parameters and plot the posterior on top of the prior.

C. Compute the posterior mean \(\mathbb{E}[\theta\mid x] = \dfrac{\alpha_{\text{post}}}{\alpha_{\text{post}}+\beta_{\text{post}}}\) and compare it to the MLE from Section 1. Which one is closer to the old 20% baseline, and why does that make sense?

# B. Posterior parameters and plot
alpha_post = alpha0 + x  # TODO: alpha0 + x
beta_post  = ta  # TODO: beta0 + (n - x)

# C. Posterior mean vs MLE
posterior_mean = None  # TODO: alpha_post / (alpha_post + beta_post)

print(f"MLE:            {x/n:.4f}")
print(f"Posterior mean: {posterior_mean:.4f}")

D. Visualize the density of the posterior distribution using beta.pdf on top of the prior distribution.

---

3. Credible Intervals and Posterior Predictive Simulation

A single number (the posterior mean) hides how uncertain we still are about true checkout rate \(\theta\). We’ll now quantify that uncertainty and use it to make a decision.

3.1. 95% Credible Interval

A. Use beta.ppf (the inverse CDF, i.e., quantile function) to find the 2.5% and 97.5% quantiles of the posterior \(\text{Beta}(\alpha_{\text{post}}, \beta_{\text{post}})\). This is your 95% credible interval: “given the data, there is a 95% (posterior) probability that \(\theta\) lies in this interval.”

# A. Analytic 95% credible interval
ci_low, ci_high = None, None  # TODO

print(f"95% Credible Interval for theta: [{ci_low:.4f}, {ci_high:.4f}]")

3.2. Confirm the Interval by Simulation

Rather than trusting the formula blindly, let’s simulate directly from the posterior and check the interval numerically — this is also good practice for the harder posterior predictive question next.

B. Draw 10,000 samples of \(\theta\) from \(\text{Beta}(\alpha_{\text{post}}, \beta_{\text{post}})\) using np.random.beta, then take the 2.5th and 97.5th percentiles (np.percentile) of your samples. Compare to part A.

# B. Simulate from the posterior
n_sim = 10_000
theta_samples = None  # TODO: simulate samples from the computed posterior distribution.

ci_low_sim = None   # TODO
ci_high_sim = None  # TODO

print(f"Simulated 95% CI: [{ci_low_sim:.4f}, {ci_high_sim:.4f}]")

3.3. Making the Business Decision

C. Using your theta_samples from part B, compute the posterior probability that the new page truly beats the old 20% baseline, i.e. \(\mathbb{P}(\theta > 0.20 \mid x)\). Hint: this is just the fraction of your samples that are above 0.20.

D. Suppose the company will get 100 new visitors next week. For each posterior sample of \(\theta\), simulate the number of purchases out of 100 using np.random.binomial(100, theta_samples). This gives you a posterior predictive distribution — it accounts for both our uncertainty about \(\theta\) and the randomness of future visitors. Plot a histogram and report the average predicted number of purchases.

E. Based on C and D, would you recommend launching the new checkout page? Justify your answer in 2-3 sentences.

# C. Posterior probability that the new page beats the baseline
prob_better = None  # TODO

print(f"P(theta > 0.20 | data) = {prob_better:.4f}")

# D. Posterior predictive simulation for 100 future visitors
future_n = 100
future_purchases = None  # TODO: for each of 10 000 samples of theta, how many would purchase among new 100 visitors?

print(f"Predicted purchases out of {future_n} visitors: mean = {future_purchases.mean():.2f}, "
      f"95% range = [{np.percentile(future_purchases, 2.5):.0f}, {np.percentile(future_purchases, 97.5):.0f}]")

plt.hist(future_purchases, bins=range(0, future_n+2), color='#d97f2e', edgecolor='white')
plt.xlabel('Number of purchases out of 100 future visitors'); plt.title('Posterior predictive distribution')
plt.show()

E. Your recommendation here:

…

---

4. Going Further

4.1. Sensitivity Analysis: Does the Prior Matter?

A natural worry with Bayesian methods is: “aren’t we just making up the prior?” Let’s check how much the choice of prior actually matters once we have data.

A. Repeat the posterior computation from Section 2 for three different priors on the same data (\(n=50\), \(x=14\)):

Prior	\(\alpha_0\)	\(\beta_0\)	Interpretation
Weak / uninformative	1	1	“I have no idea, \(\theta\) could be anything” (Uniform)
Moderate (used above)	4	16	“I believe ~20%, but not too strongly”
Strong / very confident	40	160	“I am very confident it’s ~20%”

Use the conjugate rule to get each posterior, and plot all three posteriors on one figure. What happens to the posterior mean as the prior becomes more confident? Which prior lets the data “speak the loudest”?

# A. Sensitivity to the choice of prior
priors = [(1, 1, 'Weak: Beta(1,1)'), (4, 16, 'Moderate: Beta(4,16)'), (40, 160, 'Strong: Beta(40,160)')]

# To do

B. In a markdown cell, briefly answer: if the company collected data from 5,000 visitors instead of 50 (assuming the same ~28% purchase rate), would you expect the three posteriors above to still look so different? Why or why not?

Your answer here:

…

4.2. Bonus: A Different Kind of Data — Counting Calls (Poisson-Gamma Model)

Not all data is “success/failure”. Suppose instead we track the number of calls per hour at a customer support center over 10 hours:

5, 7, 4, 6, 9, 5, 8, 6, 7, 5

We model the number of calls per hour as \(y_i \mid \lambda \sim \text{Poisson}(\lambda)\), where \(\lambda\) is the unknown average rate of calls per hour.

TipPoisson-Gamma Model

New conjugate update rule (given to you, no derivation needed):

If \(\lambda \sim \text{Gamma}(\alpha_0, \beta_0)\) (prior, rate parameterization) and we observe \(n\) hours with a total of \(S=\sum_i y_i\) calls, then:

\[(\lambda \mid y) \sim \text{Gamma}(\alpha_0 + S,\ \ \beta_0 + n) \quad \text{(posterior)}\]

The posterior mean is \(\mathbb{E}[\lambda \mid y] = \dfrac{\alpha_0+S}{\beta_0+n}\).

We’ll use a weakly informative prior \(\lambda \sim \text{Gamma}(\alpha_0=2, \beta_0=1)\) (prior mean = 2 calls/hour, easily overruled by data).

A. Compute the posterior parameters and the posterior mean rate of calls per hour. Compare it to the plain sample average of the 10 observed hours.

B. Compute a 95% credible interval for \(\lambda\) using gamma.ppf(..., a=alpha_post, scale=1/beta_post).

C. Simulate the posterior predictive distribution for the number of calls in the next hour: draw lambda samples from the posterior with np.random.gamma, then draw one Poisson call-count per sample with np.random.poisson. Plot a histogram and report \(\mathbb{P}(\text{more than 10 calls next hour})\).

# A, B, C. Poisson-Gamma model for call center data
calls = np.array([5, 7, 4, 6, 9, 5, 8, 6, 7, 5])
n_hours = len(calls)
S = calls.sum()

alpha0_g, beta0_g = 2, 1  # prior Gamma(shape=2, rate=1)

# A. Posterior parameters and mean
alpha_post_g = None  # TODO: alpha0_g + S
beta_post_g  = None  # TODO: beta0_g + n_hours
posterior_mean_lambda = None  # TODO: alpha_post_g / beta_post_g

print(f"Sample average calls/hour: {calls.mean():.4f}")
print(f"Posterior mean of lambda:  {posterior_mean_lambda:.4f}")

# B. 95% credible interval
ci_low_g, ci_high_g = None, None  # TODO: gamma.ppf([0.025, 0.975], a=alpha_post_g, scale=1/beta_post_g)
print(f"95% Credible Interval for lambda: [{ci_low_g:.4f}, {ci_high_g:.4f}]")

# C. Posterior predictive simulation for next hour
lambda_samples = None       # TODO: np.random.gamma(alpha_post_g, 1/beta_post_g, size=10_000)
calls_next_hour = None      # TODO: np.random.poisson(lambda_samples)

prob_more_than_10 = None    # TODO: np.mean(calls_next_hour > 10)
print(f"P(more than 10 calls next hour) = {prob_more_than_10:.4f}")

plt.hist(calls_next_hour, bins=range(0, 20), color='#345a8b', edgecolor='white')
plt.xlabel('Simulated calls in the next hour'); plt.title('Posterior predictive: calls center')
plt.show()

---

Wrap-up

In this lab you:

• Went from a purely data-driven MLE to a Bayesian estimate that blends prior knowledge with data (Section 1-2).
• Quantified uncertainty with a credible interval, both analytically and by simulation (Section 3.1-3.2).
• Used posterior simulation to answer a real business question and to build a posterior predictive distribution (Section 3.3).
• Saw how the strength of the prior trades off against the amount of data (Section 4.1).
• Applied the same Bayesian recipe to a completely different model — Poisson-Gamma — for count data (Section 4.2).

The recipe is always the same:

choose a model → choose a prior → update with data → summarize and simulate the posterior.