Skip to contents

Computes \(P(K_J = k \mid a, b)\) for \(k = 0, 1, \ldots, J\) when \(\alpha \sim \mathrm{Gamma}(a, b)\) (shape-rate parameterization).

Usage

pmf_K_marginal(
  J,
  a,
  b,
  logS,
  M = .QUAD_NODES_DEFAULT,
  M_verify = NULL,
  abs_tol = 1e-10,
  rel_tol = 1e-08,
  strict = FALSE
)

Arguments

J

Integer; sample size (positive integer >= 1).

a

Numeric; shape parameter of Gamma prior (> 0).

b

Numeric; rate parameter of Gamma prior (> 0).

logS

Matrix; pre-computed log-Stirling matrix from compute_log_stirling.

M

Integer; number of quadrature nodes (default: 80).

M_verify

Optional independent quadrature order satisfying the package verification rule (at least max(2*M, M+40), within the supported ceiling).

abs_tol, rel_tol

Non-negative absolute and relative tolerances for the selected-versus-verification L1 discrepancy and PMF/direct-moment audits.

strict

Logical; require successful higher-order verification when TRUE; otherwise return the selected PMF with explicit approximation metadata.

Value

Numeric vector of length \(J+1\) containing \(P(K_J = k \mid a, b)\) for \(k = 0, 1, \ldots, J\). Entry [1] corresponds to \(k=0\) and always equals 0. The vector sums to 1 and carries a "marginal_metadata" attribute describing support, normalization, quadrature, verification status, and the absence of support truncation.

Details

Uses Gauss-Laguerre quadrature to numerically evaluate: $$P(K_J = k \mid a, b) = \int_0^\infty P(K_J = k \mid \alpha) \cdot g_{a,b}(\alpha) d\alpha$$ $$\approx \sum_{m=1}^M \tilde{w}_m \cdot P(K_J = k \mid \alpha_m)$$

where \(P(K_J = k \mid \alpha)\) is the Antoniak distribution from Module 04 and \((\alpha_m, \tilde{w}_m)\) are the transformed quadrature nodes and normalized weights from Module 02.

Implementation: All mixing is performed in log-space for numerical stability. This is critical for large J or extreme parameter values. A "converged" result additionally requires higher-order PMF L1 agreement and PMF-versus-direct-moment agreement at both quadrature orders.

Key properties:

  • \(P(K_J = 0) = 0\) always (at least one cluster exists)

  • The PMF sums to 1

  • Moments from the PMF match exact_K_moments() within numerical tolerance

  • Mode is typically near \(E[K_J]\) but may differ

References

Antoniak, C. E. (1974). Mixtures of Dirichlet Processes with Applications to Bayesian Nonparametric Problems. The Annals of Statistics, 2(6), 1152-1174.

Examples

# Pre-compute Stirling numbers
logS <- compute_log_stirling(50)

# Compute marginal PMF for J=50, Gamma(1.5, 0.5) prior
pmf <- pmf_K_marginal(50, 1.5, 0.5, logS)

# Verify normalization
sum(pmf)
#> [1] 1

# Most likely number of clusters
which.max(pmf) - 1
#> [1] 6

# Compare mean with exact_K_moments
k_vals <- 0:50
mean_pmf <- sum(k_vals * pmf)
exact <- exact_K_moments(50, 1.5, 0.5)
abs(mean_pmf - exact$mean)
#> [1] 8.881784e-15