Marginal Moments of K_J under Gamma Prior

Computes the marginal mean $E[K_J]$ and variance $Var(K_J)$ when the DP concentration parameter follows a Gamma(a, b) prior.

Usage

exact_K_moments(
  J,
  a,
  b,
  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).
M: Integer; number of quadrature nodes (default: 80).
M_verify: Optional integer quadrature order at least max(2*M, M+40) and no greater than 512. When supplied, the selected-order moments are compared with a fresh higher-order computation. When omitted, the result is explicitly marked as an unverified approximation. For M > 256, no admissible verification order exists; only omitted M_verify with strict = FALSE is allowed.
abs_tol, rel_tol: Non-negative absolute and relative tolerances used for selected-versus-verification comparisons.
strict: Logical; if TRUE, require higher-order verification and raise a typed convergence error when it is absent or fails.

Value

A named list with components:

mean: Marginal mean $E[K_J | a, b]$
var: Marginal variance $Var(K_J | a, b)$
sd: Marginal standard deviation
cv: Coefficient of variation (sd/mean)
status: Either "converged" after successful higher-order verification or "approximate".
decomposition: Within-alpha and between-alpha components of the law of total variance.
quadrature: Selected/verification order, tolerances, discrepancies, and node-rule provenance.

Details

Uses Gauss-Laguerre quadrature to numerically evaluate: $$M_1(a,b) = E_{\alpha \sim \Gamma(a,b)}[\mu_J(\alpha)]$$ $$V(a,b) = E[v_J(\alpha)] + E[\mu_J(\alpha)^2] - M_1^2$$

where $\mu_J(\alpha)$ and $v_J(\alpha)$ are the conditional mean and variance from Module 03.

The Law of Total Variance decomposes the marginal variance into:

Within-alpha variance: $E[v_J(\alpha)]$
Between-alpha variance: $Var(\mu_J(\alpha))$

The function name is retained for compatibility. The returned values use a fixed Gauss–Laguerre rule and are therefore numerical approximations. Supply M_verify to obtain an explicit higher-order agreement check; the selected-order result is always returned and is never silently replaced. The quadrature metadata reports the required verification order and whether that order is available under the 512-node implementation ceiling.

Key properties:

The mean is bounded: $1 \leq E[K_J] \leq J$
Gamma mixing adds the non-negative between-$\alpha$ component to the average conditional variance; this identity alone does not imply a universal marginal variance-to-mean ratio
The marginal variance equals the average conditional variance plus the non-negative between-$\alpha$ variance of the conditional mean

References

Antoniak, C. E. (1974). Mixtures of Dirichlet Processes. The Annals of Statistics, 2(6), 1152-1174.

Examples

# Example: J=50, Gamma(1.5, 0.5) prior
result <- exact_K_moments(50, 1.5, 0.5)
print(result)
#> $mean
#> [1] 8.355487
#> 
#> $var
#> [1] 22.76895
#> 
#> $sd
#> [1] 4.771682
#> 
#> $cv
#> [1] 0.5710837
#> 
#> $status
#> [1] "approximate"
#> 
#> $decomposition
#> $decomposition$within_alpha
#> [1] 5.024824
#> 
#> $decomposition$between_alpha
#> [1] 17.74413
#> 
#> 
#> $quadrature
#> $quadrature$schema_version
#> [1] 1
#> 
#> $quadrature$engine
#> [1] "gauss-laguerre"
#> 
#> $quadrature$M_selected
#> [1] 80
#> 
#> $quadrature$M_verification
#> [1] NA
#> 
#> $quadrature$M_verification_required
#> [1] 160
#> 
#> $quadrature$verification_available
#> [1] TRUE
#> 
#> $quadrature$status
#> [1] "approximate"
#> 
#> $quadrature$reason
#> [1] "fixed_order_unverified"
#> 
#> $quadrature$verification_performed
#> [1] FALSE
#> 
#> $quadrature$verification_passed
#> [1] NA
#> 
#> $quadrature$absolute_tolerance
#> [1] 1e-10
#> 
#> $quadrature$relative_tolerance
#> [1] 1e-08
#> 
#> $quadrature$mean_difference
#> [1] NA
#> 
#> $quadrature$variance_difference
#> [1] NA
#> 
#> $quadrature$mean_tolerance
#> [1] NA
#> 
#> $quadrature$variance_tolerance
#> [1] NA
#> 
#> $quadrature$selected
#> $quadrature$selected$schema_version
#> [1] 1
#> 
#> $quadrature$selected$distribution
#> [1] "Gamma(shape, rate)"
#> 
#> $quadrature$selected$shape
#> [1] 1.5
#> 
#> $quadrature$selected$rate
#> [1] 0.5
#> 
#> $quadrature$selected$M_selected
#> [1] 80
#> 
#> $quadrature$selected$M_verification
#> [1] NA
#> 
#> $quadrature$selected$normalization
#> [1] "log-sum-exp"
#> 
#> $quadrature$selected$normalized_weight_sum
#> [1] 1
#> 
#> $quadrature$selected$node_rule
#> $quadrature$selected$node_rule$schema_version
#> [1] 1
#> 
#> $quadrature$selected$node_rule$engine
#> [1] "golub-welsch-dense-eigen"
#> 
#> $quadrature$selected$node_rule$M_selected
#> [1] 80
#> 
#> $quadrature$selected$node_rule$alpha_param
#> [1] 0.5
#> 
#> $quadrature$selected$node_rule$cache_key
#> [1] "M=80|alpha=0.5"
#> 
#> $quadrature$selected$node_rule$cache_hit
#> [1] TRUE
#> 
#> $quadrature$selected$node_rule$zero_normalized_weights
#> [1] 29
#> 
#> $quadrature$selected$node_rule$normalized_weight_sum
#> [1] 1
#> 
#> 
#> 
#> $quadrature$verification
#> NULL
#> 
#> 

# Verify the law-of-total-variance decomposition
sum(unlist(result$decomposition)) - result$var
#> [1] 0

# Verify mean bounds
1 <= result$mean && result$mean <= 50  # TRUE
#> [1] TRUE