1.1 Likelihood

Suppose we observe \(n\) independent binary trials — each a success (1) or failure (0) — with a common success probability \(\theta \in (0,1)\). The total number of successes \(y = \sum_i y_i\) follows:

\[ y \mid \theta \;\sim\; \mathrm{Binomial}(n,\, \theta). \]

The likelihood is

\[ p(y \mid \theta) \;=\; \binom{n}{y}\,\theta^y\,(1-\theta)^{n-y}, \]

and everything relevant for updating our beliefs about \(\theta\) is captured by the two sufficient statistics \(y\) (successes) and \(n - y\) (failures).

1.2 Prior

The parameter \(\theta\) is a probability, so its support is \([0,1]\). The Beta distribution is the natural conjugate prior:

\[ p(\theta) \;=\; \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a,b)}, \qquad \theta \in [0,1], \]

with shape parameters \(a > 0\) and \(b > 0\). Its mean is \(a/(a+b)\) and its variance is \(ab/[(a+b)^2(a+b+1)]\).

Interpreting the parameters. Think of \(a - 1\) as a “prior number of successes” and \(b - 1\) as a “prior number of failures” accumulated before the current study. The total “prior sample size” is \(a + b - 2\) and the prior mean is \(a/(a+b)\).

1.3 Posterior derivation

Multiply the likelihood kernel \(\theta^y(1-\theta)^{n-y}\) by the prior kernel \(\theta^{a-1}(1-\theta)^{b-1}\):

\[ p(\theta \mid y) \;\propto\; \theta^{(a+y)-1}(1-\theta)^{(b+n-y)-1}. \]

This is the kernel of a Beta distribution, so the conjugate posterior is:

\[ \boxed{\theta \mid y \;\sim\; \mathrm{Beta}(a + y,\; b + n - y).} \]

The update rule is simple:

• Add the number of successes \(y\) to the first shape parameter.
• Add the number of failures \(n - y\) to the second shape parameter.

1.4 The posterior mean as a compromise

\[ \mathbb{E}[\theta \mid y] \;=\; \frac{a + y}{a + b + n} \;=\; \frac{a+b}{a+b+n}\cdot\frac{a}{a+b} \;+\; \frac{n}{a+b+n}\cdot\frac{y}{n}. \]

This is a weighted average of the prior mean \(a/(a+b)\) and the data frequency \(y/n\), with weights proportional to the “prior sample size” \(a+b\) and the actual sample size \(n\). When \(n\) is much larger than \(a + b\) the data dominate; when \(n\) is small the prior pulls the estimate toward its own mean.

3.1 Fitting with glmb() using dBeta()

The scalar Beta–Binomial model generalizes to glmb() with family = binomial(link = "identity") and the dBeta() prior. In this intercept-only setting the single regression coefficient is \(\theta\) directly on the probability scale.

df_bech <- data.frame(y = pass)

bech_beta <- matrix(a0_b / (a0_b + b0_b), nrow = 1, ncol = 1,
                    dimnames = list(NULL, "(Intercept)"))

bech_pf <- dBeta(shape1 = a0_b, shape2 = b0_b, beta = bech_beta)

set.seed(2026)
fit_bech <- glmb(
  n       = 20000,
  y ~ 1,
  data    = df_bech,
  weights = rep(1L, n_bech),
  family  = binomial(link = "identity"),
  pfamily = bech_pf
)
print(fit_bech)
#> 
#> Call:  glmb(formula = y ~ 1, family = binomial(link = "identity"), pfamily = bech_pf, 
#>     n = 20000, data = df_bech, weights = rep(1L, n_bech))
#> 
#> Posterior Mean Coefficients:
#> (Intercept)  
#>      0.4477  
#> 
#> Effective Number of Parameters: 0.9969965 
#> Expected Residual Deviance: 2468.272 
#> DIC: 2469.269

bech_compare <- data.frame(
  Dataset = "Bechdel test",
  Posterior = bech_analytic$Posterior,
  `Analytic Mean` = bech_analytic$Mean,
  `Analytic SD`   = bech_analytic$SD,
  `glmb Post.Mean` = fit_bech$coef.means["(Intercept)"],
  `glmb Post.Sd`   = sd(fit_bech$coefficients[, "(Intercept)", drop = TRUE]),
  check.names = FALSE
)
knitr::kable(bech_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 dBeta() implements the exact conjugate update.