Nothing
Chapter 02-S04: Gamma–Poisson Conjugacy for One Count Rate"
In glmbayes: Bayesian Generalized Linear Models (IID Samples)
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)
library(glmbayes)
if (requireNamespace("bayesrules", quietly = TRUE)) library(bayesrules)
1. The model
1.1 Likelihood
Suppose we observe (n) independent counts (y_1, \dots, y_n) drawn from a Poisson distribution with unknown rate (\lambda > 0):
[
Y_i \mid \lambda \;\sim\; \mathrm{Poisson}(\lambda), \qquad i = 1,\dots,n.
]
The joint likelihood is
[
p(y_1,\dots,y_n \mid \lambda)
\;\propto\; e^{-n\lambda}\,\lambda^{\sum_i y_i}.
]
Everything relevant for updating beliefs about (\lambda) is captured by two sufficient statistics: (n) (number of observations) and (\sum_i y_i) (total observed count).
With exposure. In many applications the counts are observed over different exposure periods (e_i) (time, area, number of opportunities). The model becomes (Y_i \mid \lambda \sim \mathrm{Poisson}(e_i \lambda)), and the sufficient statistics become (\sum_i e_i) and (\sum_i y_i). An exposure-weighted example appears in Appendix A ([@Albert2009]).
1.2 Prior
The rate (\lambda > 0) calls for a positive prior. The Gamma distribution in shape–rate form is the natural conjugate choice:
[
p(\lambda) \;=\; \frac{\beta^\alpha}{\Gamma(\alpha)}\,\lambda^{\alpha-1}\,e^{-\beta\lambda},
\qquad \lambda > 0,
]
with shape (\alpha > 0) and rate (\beta > 0). Its mean is (\alpha/\beta) and its variance is (\alpha/\beta^2).
Interpreting the parameters. Think of (\alpha - 1) as a "prior total count" and (\beta) as a "prior number of observations". The prior mean (\alpha/\beta) is your pre-data estimate of the rate, and larger (\beta) encodes a stronger commitment to that estimate.
Note: many textbooks use shape–scale form ((\theta = 1/\beta)); glmbayes uses the shape–rate convention throughout.
1.3 Posterior derivation
Multiply the likelihood kernel (e^{-n\lambda}\lambda^{\sum y_i}) by the prior kernel (\lambda^{\alpha-1}e^{-\beta\lambda}):
[
p(\lambda \mid y)
\;\propto\;
\lambda^{(\alpha + \sum y_i) - 1}\,e^{-(\beta + n)\lambda}.
]
This is a Gamma kernel, so the conjugate posterior is:
[
\boxed{\lambda \mid y \;\sim\; \mathrm{Gamma}!\left(\alpha + \textstyle\sum_i y_i,\;\; \beta + n\right).}
]
Update rule: add the total count to the shape, add the sample size (or total exposure) to the rate.
The posterior mean is
[
\mathbb{E}[\lambda \mid y]
\;=\; \frac{\alpha + \sum y_i}{\beta + n}
\;=\;
\frac{\beta}{\beta + n}\cdot\frac{\alpha}{\beta}
\;+\;
\frac{n}{\beta + n}\cdot\bar{y},
]
again a weighted average of the prior mean and the data mean.
2. bayesrules illustration
The bayesrules package [@JohnsonLauEtAl2022] provides plot_gamma_poisson() and summarize_gamma_poisson() for visualizing the Gamma–Poisson update. The same scenario appears in the Bayes Rules! conjugate vignette and in Chapter 12 of the book.
Scenario. You observe daily bicycle rental counts (y_1, \dots, y_9) from the first nine days of a year. The data are (y = (1, 1, 1, 0, 0, 0, 0, 0, 0)), giving (\sum y = 3) across (n = 9) days. Your prior belief is a (\Gamma(3, 4)) distribution on the daily rate (\lambda) (prior mean = 0.75 rentals/day).
br_y_df <- data.frame(y = c(1, 1, 1, rep(0, 6))) ## n = 9, sum(y) = 3
br_n <- nrow(br_y_df)
br_shape <- 3
br_rate <- 4
## Analytic posterior
post_shape_br <- br_shape + sum(br_y_df$y) ## 3 + 3 = 6
post_rate_br <- br_rate + br_n ## 4 + 9 = 13
library(bayesrules)
plot_gamma_poisson(
shape = br_shape,
rate = br_rate,
sum_y = sum(br_y_df$y),
n = br_n,
prior = TRUE,
likelihood = TRUE,
posterior = TRUE
)
summarize_gamma_poisson(shape = br_shape, rate = br_rate,
sum_y = sum(br_y_df$y), n = br_n)
Reading the output. The posterior is (\Gamma(6, 13)) with mean (6/13 \approx 0.46). The data (mean (\bar y = 0.33)) have pulled the posterior mean down from the prior mean (0.75), and the posterior is noticeably tighter.
gp_analytic <- data.frame(
Example = "Bayes Rules! daily counts",
n = br_n,
`sum(y)` = sum(br_y_df$y),
Posterior = sprintf("Gamma(%d, %d)", post_shape_br, post_rate_br),
Mean = post_shape_br / post_rate_br,
SD = sqrt(post_shape_br) / post_rate_br,
check.names = FALSE
)
knitr::kable(gp_analytic, digits = 4,
caption = "Conjugate Gamma--Poisson posterior (Bayes Rules! scenario)")
3. Fitting with glmb()
The scalar Gamma–Poisson model maps to glmb() with family = poisson(link = "identity") and dGamma(shape, rate, beta, Inv_Dispersion = FALSE). The single coefficient is the rate (\lambda).
gp_beta <- matrix(br_shape / br_rate, 1L, 1L, dimnames = list(NULL, "(Intercept)"))
gp_pf <- dGamma(shape = br_shape, rate = br_rate,
beta = gp_beta, Inv_Dispersion = FALSE)
set.seed(2026)
fit_gp <- glmb(
n = 20000,
y ~ 1,
data = br_y_df,
weights = rep(1L, br_n),
family = poisson(link = "identity"),
pfamily = gp_pf
)
print(fit_gp)
gp_compare <- data.frame(
Example = "Bayes Rules! daily counts",
Posterior = gp_analytic$Posterior,
`Analytic Mean` = gp_analytic$Mean,
`Analytic SD` = gp_analytic$SD,
`glmb Post.Mean` = fit_gp$coef.means["(Intercept)"],
`glmb Post.Sd` = sd(fit_gp$coefficients[, "(Intercept)", drop = TRUE]),
check.names = FALSE
)
knitr::kable(gp_compare, digits = 4,
caption = "Analytic vs. glmb() posterior mean and SD")
The glmb Post.Mean and glmb Post.Sd columns should match the analytic Mean and SD to Monte Carlo error because dGamma(Inv_Dispersion = FALSE) implements the exact conjugate update.
4. Connection to glmbayes
The scalar Gamma–Poisson model is the conjugate prototype for Poisson glmb() with family = poisson(link = "identity"). In that setting the single coefficient is the rate (\lambda), the prior is dGamma(shape, rate, beta, Inv_Dispersion = FALSE), and the posterior update is the closed form above.
For general Poisson regression with covariates (log link), conjugacy is lost and glmb() uses envelope-based accept–reject sampling with a dNormal() prior on the regression coefficients (Chapter 10). The Bayes Rules! equality_index Poisson regression example is worked in Chapter 10, Appendix A.
Appendix A. Heart transplant mortality (@Albert2009)
The hearttransplants dataset [@Albert2020LearnBayes; @ChristiansenMorris1995] records, for each of 94 U.S. hospitals, the exposure e (expected deaths under the national baseline rate) and observed deaths y within 30 days of surgery. The model is (y_i \mid \lambda_i \sim \mathrm{Poisson}(e_i \lambda_i)) where (\lambda_i) is the hospital's standardized mortality ratio.
library(LearnBayes)
data("hearttransplants")
cat(sprintf("%d hospitals | sum(y) = %d | sum(e) = %.0f\n",
nrow(hearttransplants), sum(hearttransplants$y), sum(hearttransplants$e)))
Albert's prior. Albert (Ch. 3.2) uses a (\Gamma(16, 15174)) prior on (\lambda), derived from historical national data (prior mean (\approx 0.00105)). He focuses on two specific hospitals:
| Hospital | Exposure (e) | Deaths (y) | Analytic posterior | Posterior mean | Posterior SD |
|----------|-----------------|-------------|-------------------|----------------|--------------|
| A | 66 | 1 | (\Gamma(17,\; 15240)) | (17/15240 \approx 0.001115) | (\sqrt{17}/15240 \approx 0.000271) |
| B | 1767 | 4 | (\Gamma(20,\; 16941)) | (20/16941 \approx 0.001181) | (\sqrt{20}/16941 \approx 0.000264) |
alpha0_ht <- 16L; beta0_ht <- 15174L
ex_A <- 66L; yobs_A <- 1L
ex_B <- 1767L; yobs_B <- 4L
ht_albert <- data.frame(
Hospital = c("A", "B"),
e = c(ex_A, ex_B),
y = c(yobs_A, yobs_B),
Posterior = c("Gamma(17, 15240)", "Gamma(20, 16941)"),
Mean = c(
(alpha0_ht + yobs_A) / (beta0_ht + ex_A),
(alpha0_ht + yobs_B) / (beta0_ht + ex_B)
),
SD = c(
sqrt(alpha0_ht + yobs_A) / (beta0_ht + ex_A),
sqrt(alpha0_ht + yobs_B) / (beta0_ht + ex_B)
),
check.names = FALSE
)
knitr::kable(ht_albert, digits = 6,
caption = "Albert Ch. 3.2: analytic posterior mean and SD (shape--rate Gamma)")
A.1 Fitting with glmb()
The exposure enters as weights. The conjugate sampler updates (\Gamma(\alpha_0 + y,\; \beta_0 + e)) exactly.
ht_beta <- matrix(alpha0_ht / beta0_ht, 1L, 1L, dimnames = list(NULL, "(Intercept)"))
ht_pf <- dGamma(shape = alpha0_ht, rate = beta0_ht,
beta = ht_beta, Inv_Dispersion = FALSE)
set.seed(2026)
fit_A <- glmb(n = 20000, y ~ 1, data = data.frame(y = yobs_A),
weights = ex_A, family = poisson(link = "identity"), pfamily = ht_pf)
print(fit_A)
set.seed(2026)
fit_B <- glmb(n = 20000, y ~ 1, data = data.frame(y = yobs_B),
weights = ex_B, family = poisson(link = "identity"), pfamily = ht_pf)
print(fit_B)
ht_compare <- data.frame(
Hospital = c("A", "B"),
e = c(ex_A, ex_B),
y = c(yobs_A, yobs_B),
Posterior = c("Gamma(17, 15240)", "Gamma(20, 16941)"),
`Albert Mean` = ht_albert$Mean,
`Albert SD` = ht_albert$SD,
`glmb Post.Mean` = c(fit_A$coef.means["(Intercept)"],
fit_B$coef.means["(Intercept)"]),
`glmb Post.Sd` = c(sd(fit_A$coefficients[, "(Intercept)", drop = TRUE]),
sd(fit_B$coefficients[, "(Intercept)", drop = TRUE])),
check.names = FALSE
)
knitr::kable(ht_compare, digits = 6,
caption = "Albert analytic vs. glmb() posterior mean and SD")
Each fit is for one hospital: (Intercept) is the posterior draw for that hospital's rate (\lambda). The glmb Post.Mean and glmb Post.Sd columns should match Albert's Mean and SD to Monte Carlo error.
See also
- Chapter-02-S01 — introduction to Bayesian inference and conjugacy.
- Chapter-02-S02 — Normal–Normal conjugacy for one mean.
- Chapter-02-S03 — Beta–Binomial conjugacy for one proportion.
- Chapter-02-S05 — Gamma–Gamma conjugacy for a Gamma response rate.
- Chapter 10, Appendix A — Bayes Rules! Poisson regression (
equality_index).
- Chapter 18 — hierarchical Poisson mixed model (BikeSharing; block Gibbs with
rglmb).
Try the glmbayes package in your browser
Any scripts or data that you put into this service are public.
glmbayes documentation built on Aug. 5, 2026, 1:07 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( collapse = TRUE, comment = "#>" )
library(glmbayes) if (requireNamespace("bayesrules", quietly = TRUE)) library(bayesrules)
1. The model
1.1 Likelihood
Suppose we observe (n) independent counts (y_1, \dots, y_n) drawn from a Poisson distribution with unknown rate (\lambda > 0):
[ Y_i \mid \lambda \;\sim\; \mathrm{Poisson}(\lambda), \qquad i = 1,\dots,n. ]
The joint likelihood is
[ p(y_1,\dots,y_n \mid \lambda) \;\propto\; e^{-n\lambda}\,\lambda^{\sum_i y_i}. ]
Everything relevant for updating beliefs about (\lambda) is captured by two sufficient statistics: (n) (number of observations) and (\sum_i y_i) (total observed count).
With exposure. In many applications the counts are observed over different exposure periods (e_i) (time, area, number of opportunities). The model becomes (Y_i \mid \lambda \sim \mathrm{Poisson}(e_i \lambda)), and the sufficient statistics become (\sum_i e_i) and (\sum_i y_i). An exposure-weighted example appears in Appendix A ([@Albert2009]).
1.2 Prior
The rate (\lambda > 0) calls for a positive prior. The Gamma distribution in shape–rate form is the natural conjugate choice:
[ p(\lambda) \;=\; \frac{\beta^\alpha}{\Gamma(\alpha)}\,\lambda^{\alpha-1}\,e^{-\beta\lambda}, \qquad \lambda > 0, ]
with shape (\alpha > 0) and rate (\beta > 0). Its mean is (\alpha/\beta) and its variance is (\alpha/\beta^2).
Interpreting the parameters. Think of (\alpha - 1) as a "prior total count" and (\beta) as a "prior number of observations". The prior mean (\alpha/\beta) is your pre-data estimate of the rate, and larger (\beta) encodes a stronger commitment to that estimate.
Note: many textbooks use shape–scale form ((\theta = 1/\beta)); glmbayes uses the shape–rate convention throughout.
1.3 Posterior derivation
Multiply the likelihood kernel (e^{-n\lambda}\lambda^{\sum y_i}) by the prior kernel (\lambda^{\alpha-1}e^{-\beta\lambda}):
[ p(\lambda \mid y) \;\propto\; \lambda^{(\alpha + \sum y_i) - 1}\,e^{-(\beta + n)\lambda}. ]
This is a Gamma kernel, so the conjugate posterior is:
[ \boxed{\lambda \mid y \;\sim\; \mathrm{Gamma}!\left(\alpha + \textstyle\sum_i y_i,\;\; \beta + n\right).} ]
Update rule: add the total count to the shape, add the sample size (or total exposure) to the rate.
The posterior mean is
[ \mathbb{E}[\lambda \mid y] \;=\; \frac{\alpha + \sum y_i}{\beta + n} \;=\; \frac{\beta}{\beta + n}\cdot\frac{\alpha}{\beta} \;+\; \frac{n}{\beta + n}\cdot\bar{y}, ]
again a weighted average of the prior mean and the data mean.
2. bayesrules illustration
The bayesrules package [@JohnsonLauEtAl2022] provides plot_gamma_poisson() and summarize_gamma_poisson() for visualizing the Gamma–Poisson update. The same scenario appears in the Bayes Rules! conjugate vignette and in Chapter 12 of the book.
Scenario. You observe daily bicycle rental counts (y_1, \dots, y_9) from the first nine days of a year. The data are (y = (1, 1, 1, 0, 0, 0, 0, 0, 0)), giving (\sum y = 3) across (n = 9) days. Your prior belief is a (\Gamma(3, 4)) distribution on the daily rate (\lambda) (prior mean = 0.75 rentals/day).
br_y_df <- data.frame(y = c(1, 1, 1, rep(0, 6))) ## n = 9, sum(y) = 3 br_n <- nrow(br_y_df) br_shape <- 3 br_rate <- 4 ## Analytic posterior post_shape_br <- br_shape + sum(br_y_df$y) ## 3 + 3 = 6 post_rate_br <- br_rate + br_n ## 4 + 9 = 13
library(bayesrules) plot_gamma_poisson( shape = br_shape, rate = br_rate, sum_y = sum(br_y_df$y), n = br_n, prior = TRUE, likelihood = TRUE, posterior = TRUE )
summarize_gamma_poisson(shape = br_shape, rate = br_rate, sum_y = sum(br_y_df$y), n = br_n)
Reading the output. The posterior is (\Gamma(6, 13)) with mean (6/13 \approx 0.46). The data (mean (\bar y = 0.33)) have pulled the posterior mean down from the prior mean (0.75), and the posterior is noticeably tighter.
gp_analytic <- data.frame( Example = "Bayes Rules! daily counts", n = br_n, `sum(y)` = sum(br_y_df$y), Posterior = sprintf("Gamma(%d, %d)", post_shape_br, post_rate_br), Mean = post_shape_br / post_rate_br, SD = sqrt(post_shape_br) / post_rate_br, check.names = FALSE ) knitr::kable(gp_analytic, digits = 4, caption = "Conjugate Gamma--Poisson posterior (Bayes Rules! scenario)")
3. Fitting with glmb()
The scalar Gamma–Poisson model maps to glmb() with family = poisson(link = "identity") and dGamma(shape, rate, beta, Inv_Dispersion = FALSE). The single coefficient is the rate (\lambda).
gp_beta <- matrix(br_shape / br_rate, 1L, 1L, dimnames = list(NULL, "(Intercept)")) gp_pf <- dGamma(shape = br_shape, rate = br_rate, beta = gp_beta, Inv_Dispersion = FALSE) set.seed(2026) fit_gp <- glmb( n = 20000, y ~ 1, data = br_y_df, weights = rep(1L, br_n), family = poisson(link = "identity"), pfamily = gp_pf ) print(fit_gp)
gp_compare <- data.frame( Example = "Bayes Rules! daily counts", Posterior = gp_analytic$Posterior, `Analytic Mean` = gp_analytic$Mean, `Analytic SD` = gp_analytic$SD, `glmb Post.Mean` = fit_gp$coef.means["(Intercept)"], `glmb Post.Sd` = sd(fit_gp$coefficients[, "(Intercept)", drop = TRUE]), check.names = FALSE ) knitr::kable(gp_compare, digits = 4, caption = "Analytic vs. glmb() posterior mean and SD")
The glmb Post.Mean and glmb Post.Sd columns should match the analytic Mean and SD to Monte Carlo error because dGamma(Inv_Dispersion = FALSE) implements the exact conjugate update.
4. Connection to glmbayes
The scalar Gamma–Poisson model is the conjugate prototype for Poisson glmb() with family = poisson(link = "identity"). In that setting the single coefficient is the rate (\lambda), the prior is dGamma(shape, rate, beta, Inv_Dispersion = FALSE), and the posterior update is the closed form above.
For general Poisson regression with covariates (log link), conjugacy is lost and glmb() uses envelope-based accept–reject sampling with a dNormal() prior on the regression coefficients (Chapter 10). The Bayes Rules! equality_index Poisson regression example is worked in Chapter 10, Appendix A.
Appendix A. Heart transplant mortality (@Albert2009)
The hearttransplants dataset [@Albert2020LearnBayes; @ChristiansenMorris1995] records, for each of 94 U.S. hospitals, the exposure e (expected deaths under the national baseline rate) and observed deaths y within 30 days of surgery. The model is (y_i \mid \lambda_i \sim \mathrm{Poisson}(e_i \lambda_i)) where (\lambda_i) is the hospital's standardized mortality ratio.
library(LearnBayes) data("hearttransplants") cat(sprintf("%d hospitals | sum(y) = %d | sum(e) = %.0f\n", nrow(hearttransplants), sum(hearttransplants$y), sum(hearttransplants$e)))
Albert's prior. Albert (Ch. 3.2) uses a (\Gamma(16, 15174)) prior on (\lambda), derived from historical national data (prior mean (\approx 0.00105)). He focuses on two specific hospitals:
| Hospital | Exposure (e) | Deaths (y) | Analytic posterior | Posterior mean | Posterior SD | |----------|-----------------|-------------|-------------------|----------------|--------------| | A | 66 | 1 | (\Gamma(17,\; 15240)) | (17/15240 \approx 0.001115) | (\sqrt{17}/15240 \approx 0.000271) | | B | 1767 | 4 | (\Gamma(20,\; 16941)) | (20/16941 \approx 0.001181) | (\sqrt{20}/16941 \approx 0.000264) |
alpha0_ht <- 16L; beta0_ht <- 15174L ex_A <- 66L; yobs_A <- 1L ex_B <- 1767L; yobs_B <- 4L
ht_albert <- data.frame( Hospital = c("A", "B"), e = c(ex_A, ex_B), y = c(yobs_A, yobs_B), Posterior = c("Gamma(17, 15240)", "Gamma(20, 16941)"), Mean = c( (alpha0_ht + yobs_A) / (beta0_ht + ex_A), (alpha0_ht + yobs_B) / (beta0_ht + ex_B) ), SD = c( sqrt(alpha0_ht + yobs_A) / (beta0_ht + ex_A), sqrt(alpha0_ht + yobs_B) / (beta0_ht + ex_B) ), check.names = FALSE ) knitr::kable(ht_albert, digits = 6, caption = "Albert Ch. 3.2: analytic posterior mean and SD (shape--rate Gamma)")
A.1 Fitting with glmb()
The exposure enters as weights. The conjugate sampler updates (\Gamma(\alpha_0 + y,\; \beta_0 + e)) exactly.
ht_beta <- matrix(alpha0_ht / beta0_ht, 1L, 1L, dimnames = list(NULL, "(Intercept)")) ht_pf <- dGamma(shape = alpha0_ht, rate = beta0_ht, beta = ht_beta, Inv_Dispersion = FALSE) set.seed(2026) fit_A <- glmb(n = 20000, y ~ 1, data = data.frame(y = yobs_A), weights = ex_A, family = poisson(link = "identity"), pfamily = ht_pf) print(fit_A)
set.seed(2026) fit_B <- glmb(n = 20000, y ~ 1, data = data.frame(y = yobs_B), weights = ex_B, family = poisson(link = "identity"), pfamily = ht_pf) print(fit_B)
ht_compare <- data.frame( Hospital = c("A", "B"), e = c(ex_A, ex_B), y = c(yobs_A, yobs_B), Posterior = c("Gamma(17, 15240)", "Gamma(20, 16941)"), `Albert Mean` = ht_albert$Mean, `Albert SD` = ht_albert$SD, `glmb Post.Mean` = c(fit_A$coef.means["(Intercept)"], fit_B$coef.means["(Intercept)"]), `glmb Post.Sd` = c(sd(fit_A$coefficients[, "(Intercept)", drop = TRUE]), sd(fit_B$coefficients[, "(Intercept)", drop = TRUE])), check.names = FALSE ) knitr::kable(ht_compare, digits = 6, caption = "Albert analytic vs. glmb() posterior mean and SD")
Each fit is for one hospital: (Intercept) is the posterior draw for that hospital's rate (\lambda). The glmb Post.Mean and glmb Post.Sd columns should match Albert's Mean and SD to Monte Carlo error.
See also
- Chapter-02-S01 — introduction to Bayesian inference and conjugacy.
- Chapter-02-S02 — Normal–Normal conjugacy for one mean.
- Chapter-02-S03 — Beta–Binomial conjugacy for one proportion.
- Chapter-02-S05 — Gamma–Gamma conjugacy for a Gamma response rate.
- Chapter 10, Appendix A — Bayes Rules! Poisson regression (
equality_index). - Chapter 18 — hierarchical Poisson mixed model (BikeSharing; block Gibbs with
rglmb).
Try the glmbayes package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.