############################################################
# MODELING AND SIMULATION
# BAYESIAN EXERCISES 4.17, 4.18 AND 5.18
############################################################
# Install this package only once if needed:
# install.packages("bayesrules")
library(bayesrules)
## Warning: package 'bayesrules' was built under R version 4.5.3
############################################################
# EXERCISE 4.17
# DIFFERENT DATA, DIFFERENT POSTERIORS
############################################################
# Three employees start with the same prior:
# pi ~ Beta(4, 3)
#
# pi represents the probability that a user clicks on the ad.
############################################################
# 4.17 (a)
# Plot and describe the Beta(4,3) prior
############################################################
plot_beta(
alpha = 4,
beta = 3
)
# Calculate prior mean
prior_mean <- 4 / (4 + 3)
prior_mean
## [1] 0.5714286
# ----------------------------------------------------------
# EXPLANATION FOR 4.17(a)
#
# The prior distribution is:
#
# pi ~ Beta(4,3)
#
# The prior mean is:
#
# 4 / (4 + 3) = 0.571
#
# Therefore, before observing any new data, the employees
# believe that the probability of a person clicking the ad
# is approximately 57.1%.
#
# The Beta distribution is appropriate because pi is a
# probability, and probabilities must fall between 0 and 1.
#
# The prior is not extremely narrow, which means there is
# still uncertainty about the actual click probability.
# ----------------------------------------------------------
############################################################
# 4.17 (b)
# Determine each employee's posterior distribution
############################################################
# General Beta-Binomial updating rule:
#
# Prior:
# Beta(alpha, beta)
#
# Data:
# y successes out of n observations
#
# Posterior:
# Beta(alpha + y, beta + n - y)
##############################
# EMPLOYEE 1
##############################
# Employee 1:
# n = 1
# y = 0 clicks
employee1_alpha <- 4 + 0
employee1_beta <- 3 + (1 - 0)
employee1_alpha
## [1] 4
employee1_beta
## [1] 4
# Posterior:
# Beta(4,4)
##############################
# EMPLOYEE 2
##############################
# Employee 2:
# n = 10
# y = 3 clicks
employee2_alpha <- 4 + 3
employee2_beta <- 3 + (10 - 3)
employee2_alpha
## [1] 7
employee2_beta
## [1] 10
# Posterior:
# Beta(7,10)
##############################
# EMPLOYEE 3
##############################
# Employee 3:
# n = 100
# y = 20 clicks
employee3_alpha <- 4 + 20
employee3_beta <- 3 + (100 - 20)
employee3_alpha
## [1] 24
employee3_beta
## [1] 83
# Posterior:
# Beta(24,83)
# ----------------------------------------------------------
# EXPLANATION FOR 4.17(b)
#
# All three employees begin with the same Beta(4,3) prior,
# but they observe different amounts of data.
#
# Employee 1:
#
# Beta(4 + 0, 3 + 1)
# Beta(4,4)
#
# Employee 2:
#
# Beta(4 + 3, 3 + 7)
# Beta(7,10)
#
# Employee 3:
#
# Beta(4 + 20, 3 + 80)
# Beta(24,83)
#
# The posterior changes because each employee has different
# observed evidence.
#
# The more data that are collected, the more influence the
# observed evidence has on the posterior distribution.
# ----------------------------------------------------------
############################################################
# 4.17 (c)
# Plot prior, likelihood, and posterior for each employee
############################################################
##############################
# EMPLOYEE 1
##############################
plot_beta_binomial(
alpha = 4,
beta = 3,
y = 0,
n = 1
)
##############################
# EMPLOYEE 2
##############################
plot_beta_binomial(
alpha = 4,
beta = 3,
y = 3,
n = 10
)
##############################
# EMPLOYEE 3
##############################
plot_beta_binomial(
alpha = 4,
beta = 3,
y = 20,
n = 100
)
# ----------------------------------------------------------
# EXPLANATION FOR 4.17(c)
#
# The graphs display three pieces of Bayesian analysis:
#
# 1. Prior
# What we believed before seeing the new data.
#
# 2. Likelihood
# What the observed data suggest about pi.
#
# 3. Posterior
# Our updated belief after combining the prior and data.
#
# Employee 1 only studies one person.
# Because the sample size is extremely small, the data provide
# little evidence. Therefore, the posterior remains relatively
# close to the prior.
#
# Employee 2 studies 10 people, so the observed data have more
# influence. The posterior begins shifting toward the observed
# click rate of:
#
# 3 / 10 = 0.30
#
# Employee 3 studies 100 people. Their observed click rate is:
#
# 20 / 100 = 0.20
#
# Because this is a much larger sample, the data strongly
# influence the posterior. The posterior becomes concentrated
# near 0.20.
#
# This demonstrates an important Bayesian principle:
#
# More evidence usually gives the data more influence
# over the posterior.
# ----------------------------------------------------------
############################################################
# 4.17 (d)
# Summarize and compare the posterior distributions
############################################################
# Posterior means
employee1_mean <- employee1_alpha /
(employee1_alpha + employee1_beta)
employee2_mean <- employee2_alpha /
(employee2_alpha + employee2_beta)
employee3_mean <- employee3_alpha /
(employee3_alpha + employee3_beta)
employee1_mean
## [1] 0.5
employee2_mean
## [1] 0.4117647
employee3_mean
## [1] 0.2242991
# Posterior standard deviations
employee1_sd <- sqrt(
employee1_alpha * employee1_beta /
(
(employee1_alpha + employee1_beta)^2 *
(employee1_alpha + employee1_beta + 1)
)
)
employee2_sd <- sqrt(
employee2_alpha * employee2_beta /
(
(employee2_alpha + employee2_beta)^2 *
(employee2_alpha + employee2_beta + 1)
)
)
employee3_sd <- sqrt(
employee3_alpha * employee3_beta /
(
(employee3_alpha + employee3_beta)^2 *
(employee3_alpha + employee3_beta + 1)
)
)
employee1_sd
## [1] 0.1666667
employee2_sd
## [1] 0.1160016
employee3_sd
## [1] 0.04013738
# Create summary table
employee_summary <- data.frame(
Employee = c("Employee 1", "Employee 2", "Employee 3"),
Sample_Size = c(1, 10, 100),
Clicks = c(0, 3, 20),
Observed_Click_Rate = c(0/1, 3/10, 20/100),
Posterior = c("Beta(4,4)",
"Beta(7,10)",
"Beta(24,83)"),
Posterior_Mean = c(employee1_mean,
employee2_mean,
employee3_mean),
Posterior_SD = c(employee1_sd,
employee2_sd,
employee3_sd)
)
employee_summary
## Employee Sample_Size Clicks Observed_Click_Rate Posterior Posterior_Mean
## 1 Employee 1 1 0 0.0 Beta(4,4) 0.5000000
## 2 Employee 2 10 3 0.3 Beta(7,10) 0.4117647
## 3 Employee 3 100 20 0.2 Beta(24,83) 0.2242991
## Posterior_SD
## 1 0.16666667
## 2 0.11600156
## 3 0.04013738
# ----------------------------------------------------------
# EXPLANATION FOR 4.17(d)
#
# Employee 1:
#
# Posterior = Beta(4,4)
# Posterior mean = 0.50
#
# Even though the one person tested did not click, one
# observation is not enough evidence to completely change
# the prior belief.
#
#
# Employee 2:
#
# Posterior = Beta(7,10)
# Posterior mean = approximately 0.412
#
# The observed click rate was 30%, so the posterior moves
# downward from the prior mean of 57.1%.
#
#
# Employee 3:
#
# Posterior = Beta(24,83)
# Posterior mean = approximately 0.224
#
# The observed click rate was 20%. Because Employee 3 has
# 100 observations, their posterior is much closer to the
# observed data.
#
# The posterior standard deviation also becomes smaller as
# the amount of data increases.
#
# This means we become more certain about the click
# probability as more evidence is collected.
# ----------------------------------------------------------
############################################################
# EXERCISE 4.18
# A SEQUENTIAL EMPLOYEE
############################################################
# The fourth employee starts with the same prior:
#
# pi ~ Beta(4,3)
#
# Instead of collecting new data, they receive the data from
# Employees 1, 2, and 3 over three different days.
#
# Important Bayesian idea:
#
# Yesterday's posterior becomes today's prior.
#
############################################################
# 4.18 (a)
# Find posterior at the end of each day
############################################################
##############################
# ORIGINAL PRIOR
##############################
alpha0 <- 4
beta0 <- 3
##############################
# DAY 1
##############################
# Employee 1 data:
# 0 clicks out of 1 person
alpha1 <- alpha0 + 0
beta1 <- beta0 + 1
alpha1
## [1] 4
beta1
## [1] 4
# Day 1 posterior:
# Beta(4,4)
##############################
# DAY 2
##############################
# Employee 2 data:
# 3 clicks
# 7 non-clicks
alpha2 <- alpha1 + 3
beta2 <- beta1 + 7
alpha2
## [1] 7
beta2
## [1] 11
# Day 2 posterior:
# Beta(7,11)
##############################
# DAY 3
##############################
# Employee 3 data:
# 20 clicks
# 80 non-clicks
alpha3 <- alpha2 + 20
beta3 <- beta2 + 80
alpha3
## [1] 27
beta3
## [1] 91
# Day 3 posterior:
# Beta(27,91)
# Calculate means after each update
prior_mean <- alpha0 / (alpha0 + beta0)
day1_mean <- alpha1 / (alpha1 + beta1)
day2_mean <- alpha2 / (alpha2 + beta2)
day3_mean <- alpha3 / (alpha3 + beta3)
prior_mean
## [1] 0.5714286
day1_mean
## [1] 0.5
day2_mean
## [1] 0.3888889
day3_mean
## [1] 0.2288136
# ----------------------------------------------------------
# EXPLANATION FOR 4.18(a)
#
# Original prior:
#
# Beta(4,3)
#
# Mean = 4 / 7
# = approximately 0.571
#
#
# After Day 1:
#
# Beta(4,4)
#
# Mean = 4 / 8
# = 0.50
#
#
# After Day 2:
#
# Beta(7,11)
#
# Mean = 7 / 18
# = approximately 0.389
#
#
# After Day 3:
#
# Beta(27,91)
#
# Mean = 27 / 118
# = approximately 0.229
#
#
# Therefore, the employee's estimate of the click probability
# changes approximately as follows:
#
# 57.1% -> 50.0% -> 38.9% -> 22.9%
#
# As new evidence arrives, the employee continuously updates
# their knowledge.
#
# The data consistently suggest that the advertisement's
# click probability is lower than originally believed.
# ----------------------------------------------------------
############################################################
# 4.18 (b)
# Plot the prior and all three sequential posteriors
############################################################
pi_values <- seq(
from = 0,
to = 1,
length.out = 1000
)
# Calculate densities
prior_density <- dbeta(
pi_values,
shape1 = 4,
shape2 = 3
)
day1_density <- dbeta(
pi_values,
shape1 = 4,
shape2 = 4
)
day2_density <- dbeta(
pi_values,
shape1 = 7,
shape2 = 11
)
day3_density <- dbeta(
pi_values,
shape1 = 27,
shape2 = 91
)
# Plot the prior
plot(
pi_values,
prior_density,
type = "l",
lwd = 2,
ylim = c(
0,
max(
prior_density,
day1_density,
day2_density,
day3_density
)
),
xlab = "Probability of clicking (pi)",
ylab = "Density",
main = "Sequential Bayesian Updating"
)
# Add Day 1
lines(
pi_values,
day1_density,
lwd = 2,
lty = 2
)
# Add Day 2
lines(
pi_values,
day2_density,
lwd = 2,
lty = 3
)
# Add Day 3
lines(
pi_values,
day3_density,
lwd = 2,
lty = 4
)
# Add legend
legend(
"topright",
legend = c(
"Prior: Beta(4,3)",
"Day 1: Beta(4,4)",
"Day 2: Beta(7,11)",
"Day 3: Beta(27,91)"
),
lty = c(1, 2, 3, 4),
lwd = 2
)
# ----------------------------------------------------------
# EXPLANATION FOR 4.18(b)
#
# The graph shows how the employee's beliefs evolve over
# the three days.
#
# The original prior is centered around a relatively high
# click probability.
#
# After Day 1, the distribution moves slightly downward.
#
# After Day 2, it moves further downward because additional
# evidence suggests a lower click rate.
#
# After Day 3, the distribution becomes concentrated near
# approximately 23%.
#
# The distribution also becomes narrower.
#
# A narrower posterior means there is less uncertainty because
# the employee now has much more information.
#
# This is an example of sequential Bayesian learning:
#
# Prior
# +
# New evidence
# =
# Posterior
#
# Then:
#
# Posterior becomes the next prior.
# ----------------------------------------------------------
############################################################
# 4.18 (c)
# What happens if all three datasets are combined first?
############################################################
# Total clicks
total_clicks <- 0 + 3 + 20
total_clicks
## [1] 23
# Total observations
total_n <- 1 + 10 + 100
total_n
## [1] 111
# Total non-clicks
total_nonclicks <- total_n - total_clicks
total_nonclicks
## [1] 88
# Calculate posterior using all data at once
pooled_alpha <- 4 + total_clicks
pooled_beta <- 3 + total_nonclicks
pooled_alpha
## [1] 27
pooled_beta
## [1] 91
# Posterior:
# Beta(27,91)
# ----------------------------------------------------------
# EXPLANATION FOR 4.18(c)
#
# Across all three studies:
#
# Total clicks = 23
# Total observations = 111
# Total non-clicks = 88
#
# Combining the data with the Beta(4,3) prior gives:
#
# Beta(4 + 23, 3 + 88)
#
# Beta(27,91)
#
# This is exactly the same posterior obtained from sequential
# updating at the end of Day 3.
#
# Therefore, whether the employee:
#
# 1. Updates their beliefs one day at a time,
#
# OR
#
# 2. Waits and processes all the data together,
#
# the final posterior is the same.
#
# This happens because Bayesian updating uses all of the
# cumulative information contained in the data.
# ----------------------------------------------------------
############################################################
# EXERCISE 5.18
# COUNTING INSECTS
############################################################
# theta represents the insect density:
#
# average number of insects per square meter
#
# Prior information:
#
# Mean = 0.5 insects per square meter
# SD = 0.25 insects per square meter
#
# Observed data:
#
# First five square meters:
#
# 3, 2, 5, 1, 2
#
# Next 15 square meters:
#
# 0 insects in each area
#
############################################################
# 5.18 (a)
# Should we use a Normal or Poisson model?
############################################################
# Enter insect count data
insects <- c(
3, 2, 5, 1, 2,
rep(0, 15)
)
insects
## [1] 3 2 5 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
# Number of areas observed
n <- length(insects)
n
## [1] 20
# Total number of insects
sum_y <- sum(insects)
sum_y
## [1] 13
# Average insect count
sample_mean <- mean(insects)
sample_mean
## [1] 0.65
# ----------------------------------------------------------
# EXPLANATION FOR 5.18(a)
#
# A Poisson model should be used.
#
# The reason is that the observed variable is a count:
#
# 0 insects
# 1 insect
# 2 insects
# 3 insects
# etc.
#
# Counts are non-negative whole numbers.
#
# A Normal distribution would not be ideal because it is a
# continuous distribution and could theoretically produce
# impossible values such as:
#
# -1.5 insects
#
# Therefore, we model the insect counts as:
#
# Y_i | theta ~ Poisson(theta)
#
# where theta represents the underlying average number of
# insects per square meter.
# ----------------------------------------------------------
############################################################
# Find the Gamma prior parameters
############################################################
# We know:
#
# Prior mean = 0.5
# Prior SD = 0.25
prior_mean_gamma <- 0.5
prior_sd_gamma <- 0.25
# Calculate variance
prior_variance_gamma <- prior_sd_gamma^2
prior_variance_gamma
## [1] 0.0625
# For Gamma(shape, rate):
#
# Mean = shape / rate
#
# Variance = shape / rate^2
#
# Therefore:
#
# rate = mean / variance
prior_rate <- prior_mean_gamma /
prior_variance_gamma
prior_rate
## [1] 8
# shape = mean * rate
prior_shape <- prior_mean_gamma *
prior_rate
prior_shape
## [1] 4
# Therefore:
#
# theta ~ Gamma(4,8)
############################################################
# 5.18 (b)
# Plot prior, likelihood, and posterior
############################################################
# Display observed information
n
## [1] 20
sum_y
## [1] 13
sample_mean
## [1] 0.65
# Plot prior, likelihood and posterior
plot_gamma_poisson(
shape = prior_shape,
rate = prior_rate,
sum_y = sum_y,
n = n
)
# Posterior parameters
posterior_shape <- prior_shape + sum_y
posterior_rate <- prior_rate + n
posterior_shape
## [1] 17
posterior_rate
## [1] 28
# Posterior:
#
# Gamma(17,28)
# ----------------------------------------------------------
# EXPLANATION FOR 5.18(b)
#
# The prior distribution is:
#
# Gamma(4,8)
#
# Its mean is:
#
# 4 / 8 = 0.50 insects per square meter.
#
#
# The observed data contain:
#
# 13 insects across 20 square meters.
#
# Therefore, the observed sample mean is:
#
# 13 / 20 = 0.65 insects per square meter.
#
#
# Using Gamma-Poisson Bayesian updating:
#
# Posterior shape
# = 4 + 13
# = 17
#
# Posterior rate
# = 8 + 20
# = 28
#
# Therefore:
#
# theta | data ~ Gamma(17,28)
#
# The posterior lies between the prior belief of 0.50 and
# the observed sample mean of 0.65.
#
# The posterior moves toward the observed data because the
# 20 field observations provide new information.
#
# The result shows how Bayesian modeling combines:
#
# Prior knowledge
# +
# Observed evidence
# =
# Updated knowledge
# ----------------------------------------------------------
############################################################
# 5.18 (c)
# Find posterior mean and standard deviation
############################################################
# Posterior mean
posterior_mean_gamma <- posterior_shape /
posterior_rate
posterior_mean_gamma
## [1] 0.6071429
# Posterior standard deviation
#
# For Gamma(shape, rate):
#
# SD = sqrt(shape) / rate
posterior_sd_gamma <- sqrt(
posterior_shape
) / posterior_rate
posterior_sd_gamma
## [1] 0.1472538
# Compare with prior
prior_mean_gamma
## [1] 0.5
prior_sd_gamma
## [1] 0.25
posterior_mean_gamma
## [1] 0.6071429
posterior_sd_gamma
## [1] 0.1472538
# ----------------------------------------------------------
# EXPLANATION FOR 5.18(c)
#
# The posterior distribution is:
#
# Gamma(17,28)
#
# Posterior mean:
#
# 17 / 28
# = approximately 0.607
#
# Therefore, our updated estimate is approximately:
#
# 0.61 insects per square meter.
#
#
# Posterior standard deviation:
#
# sqrt(17) / 28
# = approximately 0.147
#
#
# Prior:
#
# Mean = 0.50
# SD = 0.25
#
# Posterior:
#
# Mean = approximately 0.607
# SD = approximately 0.147
#
# The posterior mean moved toward the observed field average
# of 0.65.
#
# The standard deviation decreased from 0.25 to approximately
# 0.147.
#
# This tells us that after collecting the field data, we have
# less uncertainty about the true insect density than we had
# before collecting the data.
# ----------------------------------------------------------
############################################################
# 5.18 (d)
# Interpret the posterior result for the biologist
############################################################
# Optional:
# Calculate a 95% credible interval for theta
credible_interval <- qgamma(
c(0.025, 0.975),
shape = posterior_shape,
rate = posterior_rate
)
credible_interval
## [1] 0.3536831 0.9279642
# ----------------------------------------------------------
# EXPLANATION FOR 5.18(d)
#
# After combining the biologist's previous knowledge with
# insect counts from 20 separate one-square-meter areas,
# the estimated insect density is approximately:
#
# 0.61 insects per square meter.
#
# This estimate is slightly higher than the prior expectation
# of 0.50 insects per square meter but slightly lower than
# the observed sample average of 0.65 insects per square meter.
#
# The posterior standard deviation is approximately 0.147,
# compared with the prior standard deviation of 0.25.
#
# Therefore, the field observations have reduced our
# uncertainty about the true insect density.
#
# In practical terms, the biologist can conclude that the
# insect population appears to be somewhat denser than was
# originally expected, although the overall insect density
# remains relatively low.
#
# The posterior distribution provides the biologist with an
# updated estimate that takes into account both previous
# knowledge and the newly observed field evidence.
# ----------------------------------------------------------
############################################################
# FINAL CONCEPTUAL SUMMARY
############################################################
# Exercise 4.17:
#
# Beta prior + Binomial data -> Beta posterior
#
# We are estimating a PROBABILITY.
#
# Example:
# Probability that someone clicks an advertisement.
# Exercise 4.18:
#
# Bayesian updating can occur sequentially:
#
# Prior -> Posterior -> New Prior -> New Posterior
#
# The final result is the same whether the information is
# processed sequentially or all at once.
# Exercise 5.18:
#
# Gamma prior + Poisson data -> Gamma posterior
#
# We are estimating a positive RATE or AVERAGE.
#
# Example:
# Average number of insects per square meter.
# Overall Bayesian idea:
#
# PRIOR KNOWLEDGE
# +
# OBSERVED EVIDENCE
# =
# UPDATED KNOWLEDGE
#
# The updated knowledge can then be used to support
# better-informed decisions.
############################################################