Exercise 4.17: Different Data, Different Posteriors

A shoe company developed a new internet advertisement for its latest sneaker. Three employees share the same \(\operatorname{Beta}(4,3)\) prior model for \(\pi\), where \(\pi\) is the probability that a user clicks on the advertisement.

The employees conducted three different studies:

Part 1: Prior Model

The shared prior model is:

\[ \pi\sim\operatorname{Beta}(4,3). \]

plot_beta(
  alpha = 4,
  beta = 3
)

The prior mean is:

\[ E(\pi) = \frac{4}{4+3} = \frac47 \approx0.5714. \]

The prior mode is:

\[ \operatorname{Mode}(\pi) = \frac{4-1}{4+3-2} = \frac35 = 0.60. \]

# Shared prior parameters
click_prior_alpha <- 4
click_prior_beta <- 3

# Prior mean, mode, and standard deviation
click_prior_mean <- click_prior_alpha /
  (click_prior_alpha + click_prior_beta)

click_prior_mode <- (click_prior_alpha - 1) /
  (click_prior_alpha + click_prior_beta - 2)

click_prior_sd <- sqrt(
  (click_prior_alpha * click_prior_beta) /
    (
      (click_prior_alpha + click_prior_beta)^2 *
      (click_prior_alpha + click_prior_beta + 1)
    )
)

click_prior_mean
## [1] 0.5714286
click_prior_mode
## [1] 0.6
click_prior_sd
## [1] 0.1749636

The prior mean is approximately \(0.5714\), the mode is \(0.60\), and the standard deviation is approximately \(0.1750\).

Before observing the data, the employees believed that the probability of a user clicking the advertisement was most likely around 60%. The prior allows uncertainty around this value but gives more support to click probabilities above 50% than to very low click probabilities.

Part 2: Posterior Models Constructed from Scratch

The shared Beta prior density has the following proportional form:

\[ f(\pi) \propto \pi^{4-1}(1-\pi)^{3-1} = \pi^3(1-\pi)^2. \]

For binomial data with \(y\) clicks among \(n\) users, the likelihood is:

\[ L(\pi\mid y) \propto \pi^y(1-\pi)^{n-y}. \]

The posterior is proportional to the prior multiplied by the likelihood:

\[ f(\pi\mid y) \propto f(\pi)L(\pi\mid y). \]

Employee 1

Employee 1 observed:

\[ n_1=1,\qquad y_1=0. \]

The likelihood is:

\[ L_1(\pi)\propto(1-\pi). \]

Multiplying the prior by the likelihood:

\[ f(\pi\mid y_1) \propto \pi^3(1-\pi)^2(1-\pi) = \pi^3(1-\pi)^3. \]

Therefore:

\[ \boxed{\pi\mid y_1\sim\operatorname{Beta}(4,4)}. \]

Employee 2

Employee 2 observed:

\[ n_2=10,\qquad y_2=3. \]

The likelihood is:

\[ L_2(\pi) \propto \pi^3(1-\pi)^7. \]

Multiplying the prior by the likelihood:

\[ f(\pi\mid y_2) \propto \pi^3(1-\pi)^2 \pi^3(1-\pi)^7 = \pi^6(1-\pi)^9. \]

Therefore:

\[ \boxed{\pi\mid y_2\sim\operatorname{Beta}(7,10)}. \]

Employee 3

Employee 3 observed:

\[ n_3=100,\qquad y_3=20. \]

The likelihood is:

\[ L_3(\pi) \propto \pi^{20}(1-\pi)^{80}. \]

Multiplying the prior by the likelihood:

\[ f(\pi\mid y_3) \propto \pi^3(1-\pi)^2 \pi^{20}(1-\pi)^{80} = \pi^{23}(1-\pi)^{82}. \]

Therefore:

\[ \boxed{\pi\mid y_3\sim\operatorname{Beta}(24,83)}. \]

Part 3: Prior, Likelihood, and Posterior Plots

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
)

The plots show that the influence of the likelihood increases as the sample size increases. Employee 1’s single observation provides little information, while Employee 3’s 100 observations provide strong evidence about \(\pi\).

Part 4: Posterior Summaries and Comparison

# Required Beta-Binomial summaries
employee_1_summary <- summarize_beta_binomial(
  alpha = 4,
  beta = 3,
  y = 0,
  n = 1
)

employee_2_summary <- summarize_beta_binomial(
  alpha = 4,
  beta = 3,
  y = 3,
  n = 10
)

employee_3_summary <- summarize_beta_binomial(
  alpha = 4,
  beta = 3,
  y = 20,
  n = 100
)

employee_1_summary
##       model alpha beta      mean mode        var        sd
## 1     prior     4    3 0.5714286  0.6 0.03061224 0.1749636
## 2 posterior     4    4 0.5000000  0.5 0.02777778 0.1666667
employee_2_summary
##       model alpha beta      mean mode        var        sd
## 1     prior     4    3 0.5714286  0.6 0.03061224 0.1749636
## 2 posterior     7   10 0.4117647  0.4 0.01345636 0.1160016
employee_3_summary
##       model alpha beta      mean      mode         var         sd
## 1     prior     4    3 0.5714286 0.6000000 0.030612245 0.17496355
## 2 posterior    24   83 0.2242991 0.2190476 0.001611009 0.04013738
# Posterior parameters
click_posterior_alpha <- c(4, 7, 24)
click_posterior_beta <- c(4, 10, 83)

# Calculate posterior statistics
click_posterior_means <- click_posterior_alpha /
  (click_posterior_alpha + click_posterior_beta)

click_posterior_modes <- (click_posterior_alpha - 1) /
  (click_posterior_alpha + click_posterior_beta - 2)

click_posterior_sds <- sqrt(
  (click_posterior_alpha * click_posterior_beta) /
    (
      (click_posterior_alpha + click_posterior_beta)^2 *
      (click_posterior_alpha + click_posterior_beta + 1)
    )
)

# Create comparison table
click_comparison <- data.frame(
  Employee = c("Employee 1", "Employee 2", "Employee 3"),
  Sample_Size = c(1, 10, 100),
  Clicks = c(0, 3, 20),
  Sample_Proportion = c(0, 0.30, 0.20),
  Posterior_Model = c(
    "Beta(4, 4)",
    "Beta(7, 10)",
    "Beta(24, 83)"
  ),
  Posterior_Mean = round(click_posterior_means, 4),
  Posterior_Mode = round(click_posterior_modes, 4),
  Posterior_SD = round(click_posterior_sds, 4)
)

knitr::kable(
  click_comparison,
  caption = "Comparison of the Employees' Posterior Models"
)
Comparison of the Employees’ Posterior Models
Employee Sample_Size Clicks Sample_Proportion Posterior_Model Posterior_Mean Posterior_Mode Posterior_SD
Employee 1 1 0 0.0 Beta(4, 4) 0.5000 0.500 0.1667
Employee 2 10 3 0.3 Beta(7, 10) 0.4118 0.400 0.1160
Employee 3 100 20 0.2 Beta(24, 83) 0.2243 0.219 0.0401

The posterior results are:

Employee Data Sample proportion Posterior Mean Mode SD
Employee 1 0 of 1 0.00 Beta(4,4) 0.5000 0.5000 0.1667
Employee 2 3 of 10 0.30 Beta(7,10) 0.4118 0.4000 0.1160
Employee 3 20 of 100 0.20 Beta(24,83) 0.2243 0.2190 0.0401

All three employees started with the same prior, but their posterior models differ because they collected different amounts of data and observed different click rates.

Employee 1 collected only one observation. This provides very little evidence, so the posterior remains relatively close to the prior and has the greatest uncertainty.

Employee 2 observed a click rate of 30% among 10 people. The posterior moves toward 30%, but the prior still has a noticeable influence because the sample is relatively small.

Employee 3 observed a click rate of 20% among 100 people. This sample provides much more information than the prior. As a result, the posterior mean is close to 20%, and the posterior standard deviation is the smallest.

Therefore, a larger sample causes the posterior model to more closely reflect the observed data and reduces uncertainty about the true click probability.

Exercise 5.18: Counting Insects

A biologist wants to estimate the density of a certain insect in the local region. Let \(\theta\) represent the expected number of insects per square meter.

The biologist’s prior understanding is represented by a Gamma model with:

\[ E(\theta)=0.5 \]

and

\[ SD(\theta)=0.25. \]

Twenty separate one-square-meter areas were inspected. The observed insect counts were:

\[ 3,\ 2,\ 5,\ 1,\ 2,\ 0,\ldots,\ 0, \]

where the final 15 areas contained no insects.

# Insect counts in 20 one-square-meter areas
insect_counts <- c(3, 2, 5, 1, 2, rep(0, 15))

# Summarize the observed data
n_areas <- length(insect_counts)
total_insects <- sum(insect_counts)
insect_sample_mean <- mean(insect_counts)
insect_sample_variance <- var(insect_counts)

n_areas
## [1] 20
total_insects
## [1] 13
insect_sample_mean
## [1] 0.65
insect_sample_variance
## [1] 1.818421

The survey recorded 13 insects across 20 square-meter areas. The observed mean was:

\[ \bar y = \frac{13}{20} = 0.65 \]

insects per square meter.

Part 1: Selecting a Sampling Model

A Poisson model is more appropriate than a Normal model.

The response variable is a count of insects within a fixed unit of area. The possible observations are nonnegative integers:

\[ 0,1,2,3,\ldots \]

A Normal model is continuous and could assign probability to impossible negative insect counts. A Poisson model is designed for counts occurring within a fixed amount of space or time.

The sampling model is:

\[ Y_i\mid\theta \sim \operatorname{Poisson}(\theta), \]

where \(Y_i\) represents the number of insects in the \(i\)-th square-meter area and \(\theta\) is the expected number of insects per square meter.

This model assumes that the areas are observed independently, are equal in size, and have approximately the same underlying insect density.

Part 2: Prior, Likelihood, and Posterior Models

Prior Gamma Model

We use a Gamma distribution with shape parameter \(\alpha\) and rate parameter \(\beta\):

\[ \theta\sim\operatorname{Gamma}(\alpha,\beta). \]

Under the shape-rate parameterization:

\[ E(\theta)=\frac{\alpha}{\beta} \]

and

\[ SD(\theta)=\frac{\sqrt{\alpha}}{\beta}. \]

The parameters can be calculated from the prior mean \(\mu\) and standard deviation \(\sigma\):

\[ \alpha = \left(\frac{\mu}{\sigma}\right)^2 \]

and

\[ \beta = \frac{\mu}{\sigma^2}. \]

Using \(\mu=0.5\) and \(\sigma=0.25\):

\[ \alpha = \left(\frac{0.5}{0.25}\right)^2 = 4 \]

and

\[ \beta = \frac{0.5}{0.25^2} = 8. \]

Therefore, the prior model is:

\[ \boxed{\theta\sim\operatorname{Gamma}(4,8)}. \]

# Prior mean and standard deviation
insect_prior_mean_given <- 0.5
insect_prior_sd_given <- 0.25

# Convert prior information into Gamma shape and rate
insect_prior_alpha <-
  (insect_prior_mean_given / insect_prior_sd_given)^2

insect_prior_beta <-
  insect_prior_mean_given / insect_prior_sd_given^2

insect_prior_alpha
## [1] 4
insect_prior_beta
## [1] 8
# Verify the prior mean and standard deviation
insect_prior_alpha / insect_prior_beta
## [1] 0.5
sqrt(insect_prior_alpha) / insect_prior_beta
## [1] 0.25

Likelihood Function

For independent Poisson observations, the likelihood is:

\[ L(\theta\mid\mathbf y) \propto \theta^{\sum y_i}e^{-n\theta}. \]

Because \(\sum y_i=13\) and \(n=20\):

\[ L(\theta\mid\mathbf y) \propto \theta^{13}e^{-20\theta}. \]

Posterior Model

The Gamma prior density is proportional to:

\[ f(\theta) \propto \theta^{\alpha-1}e^{-\beta\theta}. \]

For the \(\operatorname{Gamma}(4,8)\) prior:

\[ f(\theta) \propto \theta^3e^{-8\theta}. \]

The posterior density is proportional to the prior multiplied by the likelihood:

\[ f(\theta\mid\mathbf y) \propto \left(\theta^3e^{-8\theta}\right) \left(\theta^{13}e^{-20\theta}\right). \]

Therefore:

\[ f(\theta\mid\mathbf y) \propto \theta^{16}e^{-28\theta}. \]

This is the kernel of a Gamma distribution with shape 17 and rate 28:

\[ \boxed{\theta\mid\mathbf y\sim\operatorname{Gamma}(17,28)}. \]

# Calculate the posterior parameters
insect_posterior_alpha <-
  insect_prior_alpha + total_insects

insect_posterior_beta <-
  insect_prior_beta + n_areas

insect_posterior_alpha
## [1] 17
insect_posterior_beta
## [1] 28

Prior, Likelihood, and Posterior Plot

# Possible values of insect density
theta_values <- seq(0, 1.6, length.out = 1000)

# Prior density
insect_prior_density <- dgamma(
  theta_values,
  shape = insect_prior_alpha,
  rate = insect_prior_beta
)

# Calculate log-likelihood
insect_log_likelihood <-
  total_insects *
  log(pmax(theta_values, .Machine$double.eps)) -
  n_areas * theta_values

# Rescale the likelihood for comparison
insect_likelihood_scaled <- exp(
  insect_log_likelihood - max(insect_log_likelihood)
)

# Posterior density
insect_posterior_density <- dgamma(
  theta_values,
  shape = insect_posterior_alpha,
  rate = insect_posterior_beta
)

# Place likelihood on a comparable vertical scale
insect_likelihood_for_plot <-
  insect_likelihood_scaled *
  max(insect_prior_density, insect_posterior_density)

# Plot prior
plot(
  theta_values,
  insect_prior_density,
  type = "l",
  lwd = 3,
  col = "steelblue",
  ylim = c(
    0,
    max(
      insect_prior_density,
      insect_likelihood_for_plot,
      insect_posterior_density
    )
  ),
  main = "Prior, Likelihood, and Posterior for Insect Density",
  xlab = expression(
    theta~"(insects per square meter)"
  ),
  ylab = "Density / Rescaled likelihood"
)

# Add likelihood
lines(
  theta_values,
  insect_likelihood_for_plot,
  lwd = 3,
  col = "darkgreen",
  lty = 2
)

# Add posterior
lines(
  theta_values,
  insect_posterior_density,
  lwd = 3,
  col = "orange"
)

# Mark prior and observed means
abline(
  v = insect_prior_mean_given,
  col = "steelblue",
  lty = 3
)

abline(
  v = insect_sample_mean,
  col = "darkgreen",
  lty = 3
)

legend(
  "topright",
  legend = c(
    "Prior: Gamma(4, 8)",
    "Rescaled likelihood",
    "Posterior: Gamma(17, 28)",
    "Prior mean = 0.50",
    "Observed mean = 0.65"
  ),
  col = c(
    "steelblue",
    "darkgreen",
    "orange",
    "steelblue",
    "darkgreen"
  ),
  lty = c(1, 2, 1, 3, 3),
  lwd = c(3, 3, 3, 1, 1),
  cex = 0.8
)

The prior distribution is centered at 0.50 insects per square meter, while the observed data have a mean of 0.65. The posterior distribution lies between these two values but is closer to the observed data.

The prior rate parameter of 8 can be interpreted as approximately eight square meters of prior exposure. The survey provides 20 square meters of new exposure. Therefore, the field data provide more information and have a stronger influence on the posterior.

The posterior distribution is also narrower than the prior, indicating that the new observations reduced uncertainty about \(\theta\).

Part 3: Posterior Mean and Standard Deviation

For a Gamma distribution with shape \(\alpha\) and rate \(\beta\):

\[ E(\theta)=\frac{\alpha}{\beta} \]

and

\[ SD(\theta)=\frac{\sqrt{\alpha}}{\beta}. \]

For the posterior \(\operatorname{Gamma}(17,28)\):

\[ E(\theta\mid\mathbf y) = \frac{17}{28} \approx0.6071. \]

The posterior standard deviation is:

\[ SD(\theta\mid\mathbf y) = \frac{\sqrt{17}}{28} \approx0.1473. \]

# Posterior mean
insect_posterior_mean <-
  insect_posterior_alpha / insect_posterior_beta

# Posterior standard deviation
insect_posterior_sd <-
  sqrt(insect_posterior_alpha) /
  insect_posterior_beta

insect_posterior_mean
## [1] 0.6071429
insect_posterior_sd
## [1] 0.1472538
# Create summary table
insect_summary <- data.frame(
  Measurement = c(
    "Posterior model",
    "Posterior mean",
    "Posterior standard deviation"
  ),
  Result = c(
    "Gamma(17, 28)",
    round(insect_posterior_mean, 4),
    round(insect_posterior_sd, 4)
  )
)

knitr::kable(
  insect_summary,
  caption = "Posterior Summary for Insect Density"
)
Posterior Summary for Insect Density
Measurement Result
Posterior model Gamma(17, 28)
Posterior mean 0.6071
Posterior standard deviation 0.1473

Therefore:

  • Posterior model: \(\operatorname{Gamma}(17,28)\)
  • Posterior mean: approximately 0.6071 insects per square meter
  • Posterior standard deviation: approximately 0.1473 insects per square meter

Part 4: Conclusions for the Biologist

After combining the biologist’s prior information with observations from 20 square-meter areas, the estimated average insect density is approximately:

\[ \boxed{0.6071\text{ insects per square meter}}. \]

The posterior standard deviation is approximately:

\[ \boxed{0.1473\text{ insects per square meter}}. \]

The posterior estimate is higher than the prior mean of 0.50 but slightly lower than the observed sample mean of 0.65 insects per square meter. It represents a compromise between the prior knowledge and the field data, with the field data having the greater influence.

The insect appears to occur at a relatively low average density. However, all 13 insects were found in only the first five areas, while the remaining 15 areas contained no insects. This pattern suggests that the insects may not be evenly distributed throughout the field. They may be concentrated in particular locations or habitat patches.

The Gamma-Poisson analysis estimates the overall average density, but the biologist should consider collecting more observations across different locations. Habitat characteristics such as vegetation, moisture, soil conditions, and available food could also be recorded. This would help determine whether the numerous zero counts result from ordinary Poisson variation or from spatial clustering that requires a more flexible model.