The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Chapter 02-S04: Gamma–Poisson Conjugacy for One Count Rate

Kjell Nygren

2026-06-21

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 ((Albert 2009)).

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 (Johnson et al. 2022) 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)
#>       model shape rate      mean      mode        var        sd
#> 1     prior     3    4 0.7500000 0.5000000 0.18750000 0.4330127
#> 2 posterior     6   13 0.4615385 0.3846154 0.03550296 0.1884223

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)")

Conjugate Gamma–Poisson posterior (Bayes Rules! scenario)
Example	n	sum(y)	Posterior	Mean	SD
Bayes Rules! daily counts	9	3	Gamma(6, 13)	0.4615	0.1884

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)
#> 
#> Call:  glmb(formula = y ~ 1, family = poisson(link = "identity"), pfamily = gp_pf, 
#>     n = 20000, data = br_y_df, weights = rep(1L, br_n))
#> 
#> Posterior Mean Coefficients:
#> (Intercept)  
#>      0.4641  
#> 
#> Effective Number of Parameters: 0.5025621 
#> Expected Residual Deviance: 7.462233 
#> DIC: 13.96479

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")

Analytic vs. glmb() posterior mean and SD
	Example	Posterior	Analytic Mean	Analytic SD	glmb Post.Mean	glmb Post.Sd
(Intercept)	Bayes Rules! daily counts	Gamma(6, 13)	0.4615	0.1884	0.4641	0.1875

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 (Albert (2009))

The hearttransplants dataset (Albert 2020; Christiansen and Morris 1995) 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)))
#> 94 hospitals  |  sum(y) = 277  |  sum(e) = 294681

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)")

Albert Ch. 3.2: analytic posterior mean and SD (shape–rate Gamma)
Hospital	e	y	Posterior	Mean	SD
A	66	1	Gamma(17, 15240)	0.001115	0.000271
B	1767	4	Gamma(20, 16941)	0.001181	0.000264

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)
#> 
#> Call:  glmb(formula = y ~ 1, family = poisson(link = "identity"), pfamily = ht_pf, 
#>     n = 20000, data = data.frame(y = yobs_A), weights = ex_A)
#> 
#> Posterior Mean Coefficients:
#> (Intercept)  
#>    0.001121  
#> 
#> Effective Number of Parameters: 3.864282 
#> Expected Residual Deviance: 768.8068 
#> DIC: 904.6711

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)
#> 
#> Call:  glmb(formula = y ~ 1, family = poisson(link = "identity"), pfamily = ht_pf, 
#>     n = 20000, data = data.frame(y = yobs_B), weights = ex_B)
#> 
#> Posterior Mean Coefficients:
#> (Intercept)  
#>    0.001185  
#> 
#> Effective Number of Parameters: 351.3413 
#> Expected Residual Deviance: 101062.1 
#> DIC: 107184

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")

Albert analytic vs. glmb() posterior mean and SD
Hospital	e	y	Posterior	Albert Mean	Albert SD	glmb Post.Mean	glmb Post.Sd
A	66	1	Gamma(17, 15240)	0.001115	0.000271	0.001121	0.000270
B	1767	4	Gamma(20, 16941)	0.001181	0.000264	0.001185	0.000263

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.

References

Albert, Jim. 2009. Bayesian Computation with R. 2nd ed. Use r! Springer.

Albert, Jim. 2020. LearnBayes: Functions for Learning Bayesian Inference. https://CRAN.R-project.org/package=LearnBayes.

Christiansen, Cindy L., and Carl N. Morris. 1995. Fitting and Checking a Two-Level Poisson Model: Modeling Patient Mortality Rates in Heart Transplant Patients. (New York).

Johnson, Alicia A., Miles Q. Ott, and Mine Dogucu. 2022. Bayes Rules! An Introduction to Applied Bayesian Modeling. CRC Press. https://www.bayesrulesbook.com.

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.