Skip to contents

Compute Marginal Moments and Jacobian Simultaneously

moments_with_jacobian.Rd

Computes quadrature approximations to the marginal moments \(M_1 = E[K_J]\) and \(V = Var(K_J)\) along with the Jacobian matrix of the moment map \(F(a,b) = (M_1, V)\) using score function identities.

Usage

moments_with_jacobian(J, a, b, M = .QUAD_NODES_DEFAULT)

Arguments

J

Integer; sample size (number of observations/sites).

a

Numeric; shape parameter of the Gamma prior on \(\alpha\) (> 0).

b

Numeric; rate parameter of the Gamma prior on \(\alpha\) (> 0).

M

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

Value

A named list with components:

mean

Marginal mean \(E[K_J]\)

var

Marginal variance \(Var(K_J)\)

jacobian

2x2 Jacobian matrix with structure: $$J_F = \begin{bmatrix} \partial M_1/\partial a & \partial M_1/\partial b \\ \partial V/\partial a & \partial V/\partial b \end{bmatrix}$$

derivative_diagnostics

Order-refinement status, tolerances, and quadrature orders used for the score derivatives.

conditioning

Scale-aware singularity diagnostics for both shape-rate and solver log-parameterizations.

Details

This function uses the score identity (Lee, 2026, Section 3.2, Corollary 1) to compute derivatives without finite differences: $$\frac{\partial}{\partial\theta} E[f(\alpha)] = E[f(\alpha) \cdot s_\theta(\alpha)]$$

The Jacobian components are computed as:

  • \(\partial M_1/\partial \theta = E[(\mu_J(\alpha)-1) \cdot s_\theta(\alpha)]\)

  • \(\partial V/\partial \theta = \partial E[v_J]/\partial \theta + \partial E[\mu_J^2]/\partial \theta - 2 M_1 \partial M_1/\partial \theta\)

Subtracting the exact zero-concentration limits is an exact control variate because \(E[s_\theta(\alpha)]=0\). It removes the spurious constant times log(alpha) term that otherwise dominates finite-order shape-score quadrature. Derivatives are recomputed at a distinct higher order when the supported ceiling permits; the returned diagnostic is converged only when a full scheduled refinement is available and the two orders meet the release mixed error budget 5e-5 + 2e-3 * max(abs(selected), abs(verification)) componentwise.

Numerical Considerations:

  • The score function s_a contains log(alpha), which causes slower quadrature convergence compared to moment computation.

  • The score calculation uses at least 80 nodes and a distinct higher order when available; both orders are retained in diagnostics.

  • Conditional moments use stable finite sums without a small-alpha threshold or negative clipping.

  • The returned marginal variance uses the non-negative within/between decomposition rather than cancellation-prone second-moment subtraction.

  • Conditioning status is based on the solver-relevant log-parameter Jacobian, not on an unscaled determinant alone.

References

Lee, J. (2026). Design-Conditional Prior Elicitation for Dirichlet Process Mixtures. arXiv preprint arXiv:2602.06301.

See also

exact_K_moments for moments only, score_a, score_b for score functions

Examples

# Compute moments and Jacobian for J=50, a=2, b=1
result <- moments_with_jacobian(J = 50, a = 2.0, b = 1.0)
print(result$mean)      # E[K_J]
#> [1] 6.639693
print(result$var)       # Var(K_J)
#> [1] 12.9545
print(result$jacobian)  # 2x2 Jacobian matrix
#>           da         db
#> dM1 2.245520  -4.135585
#> dV  2.944639 -13.038228

# Use in Newton iteration
target <- c(5.0, 8.0)  # Target (E[K], Var(K))
current <- c(result$mean, result$var)
residual <- current - target
delta <- solve(result$jacobian, -residual)