Skip to contents

A1 Approximation Moment Errors

a1_moment_error.Rd

Computes the discrepancy between the A1 (shifted NegBin) approximation and exact marginal moments of \(K_J\).

Usage

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

Arguments

J

Integer; sample size.

a, b

Numeric; Gamma hyperparameters (shape, rate).

cJ

Numeric; scaling constant (default: log(J)).

M

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

M_verify

Optional independent quadrature order for the exact moments. It must satisfy the package-wide verification-order contract.

abs_tol, rel_tol

Non-negative selected-versus-verification tolerances.

strict

Logical; require verified exact-moment quadrature.

Value

A list with components:

exact_mean, exact_var

Exact moments via Gauss-Laguerre quadrature

a1_mean, a1_var

A1 (shifted NegBin) approximation moments

error_mean_abs, error_var_abs

Absolute errors

error_mean_rel, error_var_rel

Relative errors (percentage)

Details

The A1 approximation models \(K_J \approx 1 + \text{NegBin}(a, p_J)\) where \(p_J = b / (b + c_J)\).

The NegBin(a, p) moments are:

  • Mean: \(a(1-p)/p\)

  • Variance: \(a(1-p)/p^2\)

See also

exact_K_moments, DPprior_error_bounds

Examples

# Moment errors for J=50, Gamma(2, 1) prior
errors <- a1_moment_error(J = 50, a = 2, b = 1)
print(errors)
#> $exact_mean
#> [1] 6.639693
#> 
#> $exact_var
#> [1] 12.9545
#> 
#> $a1_mean
#> [1] 8.824046
#> 
#> $a1_var
#> [1] 38.43189
#> 
#> $error_mean_abs
#> [1] 2.184353
#> 
#> $error_var_abs
#> [1] 25.47739
#> 
#> $error_mean_rel
#> [1] 32.89841
#> 
#> $error_var_rel
#> [1] 196.6682
#> 
#> attr(,"a1_error_metadata")
#> attr(,"a1_error_metadata")$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$status
#> [1] "approximate"
#> 
#> attr(,"a1_error_metadata")$verified
#> [1] FALSE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature
#> attr(,"a1_error_metadata")$exact_moment_quadrature$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$engine
#> [1] "gauss-laguerre"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_verification
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_verification_required
#> [1] 160
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_available
#> [1] TRUE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$status
#> [1] "approximate"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$reason
#> [1] "fixed_order_unverified"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_performed
#> [1] FALSE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_passed
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$absolute_tolerance
#> [1] 1e-10
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$relative_tolerance
#> [1] 1e-08
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$mean_difference
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$variance_difference
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$mean_tolerance
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$variance_tolerance
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$distribution
#> [1] "Gamma(shape, rate)"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$shape
#> [1] 2
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$rate
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$M_verification
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$normalization
#> [1] "log-sum-exp"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$normalized_weight_sum
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$engine
#> [1] "golub-welsch-dense-eigen"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$alpha_param
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$cache_key
#> [1] "M=80|alpha=1"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$cache_hit
#> [1] TRUE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$zero_normalized_weights
#> [1] 29
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$normalized_weight_sum
#> [1] 1
#> 
#> 
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification
#> NULL
#> 
#> 
#> attr(,"a1_error_metadata")$proxy
#> attr(,"a1_error_metadata")$proxy$conditional
#> [1] "Poisson(alpha * cJ)"
#> 
#> attr(,"a1_error_metadata")$proxy$marginal
#> [1] "1 + NegativeBinomial(a, b/(b+cJ))"
#> 
#> attr(,"a1_error_metadata")$proxy$cJ
#> [1] 3.912023
#> 
#> 

# Compare A1 accuracy at different J values
sapply(c(25, 50, 100, 200), function(J) {
  err <- a1_moment_error(J, a = 2, b = 1)
  c(mean_err = err$error_mean_rel, var_err = err$error_var_rel)
})
#>               [,1]      [,2]      [,3]      [,4]
#> mean_err  39.18536  32.89841  27.97236  24.15421
#> var_err  263.32276 196.66824 152.98144 123.40937

# Compare different scaling constants
a1_moment_error(J = 50, a = 2, b = 1, cJ = log(50))
#> $exact_mean
#> [1] 6.639693
#> 
#> $exact_var
#> [1] 12.9545
#> 
#> $a1_mean
#> [1] 8.824046
#> 
#> $a1_var
#> [1] 38.43189
#> 
#> $error_mean_abs
#> [1] 2.184353
#> 
#> $error_var_abs
#> [1] 25.47739
#> 
#> $error_mean_rel
#> [1] 32.89841
#> 
#> $error_var_rel
#> [1] 196.6682
#> 
#> attr(,"a1_error_metadata")
#> attr(,"a1_error_metadata")$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$status
#> [1] "approximate"
#> 
#> attr(,"a1_error_metadata")$verified
#> [1] FALSE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature
#> attr(,"a1_error_metadata")$exact_moment_quadrature$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$engine
#> [1] "gauss-laguerre"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_verification
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_verification_required
#> [1] 160
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_available
#> [1] TRUE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$status
#> [1] "approximate"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$reason
#> [1] "fixed_order_unverified"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_performed
#> [1] FALSE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_passed
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$absolute_tolerance
#> [1] 1e-10
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$relative_tolerance
#> [1] 1e-08
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$mean_difference
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$variance_difference
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$mean_tolerance
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$variance_tolerance
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$distribution
#> [1] "Gamma(shape, rate)"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$shape
#> [1] 2
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$rate
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$M_verification
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$normalization
#> [1] "log-sum-exp"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$normalized_weight_sum
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$engine
#> [1] "golub-welsch-dense-eigen"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$alpha_param
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$cache_key
#> [1] "M=80|alpha=1"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$cache_hit
#> [1] TRUE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$zero_normalized_weights
#> [1] 29
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$normalized_weight_sum
#> [1] 1
#> 
#> 
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification
#> NULL
#> 
#> 
#> attr(,"a1_error_metadata")$proxy
#> attr(,"a1_error_metadata")$proxy$conditional
#> [1] "Poisson(alpha * cJ)"
#> 
#> attr(,"a1_error_metadata")$proxy$marginal
#> [1] "1 + NegativeBinomial(a, b/(b+cJ))"
#> 
#> attr(,"a1_error_metadata")$proxy$cJ
#> [1] 3.912023
#> 
#> 
a1_moment_error(J = 50, a = 2, b = 1, cJ = digamma(50) + 0.5772)
#> $exact_mean
#> [1] 6.639693
#> 
#> $exact_var
#> [1] 12.9545
#> 
#> $a1_mean
#> [1] 9.958379
#> 
#> $a1_var
#> [1] 49.08466
#> 
#> $error_mean_abs
#> [1] 3.318686
#> 
#> $error_var_abs
#> [1] 36.13016
#> 
#> $error_mean_rel
#> [1] 49.98253
#> 
#> $error_var_rel
#> [1] 278.9004
#> 
#> attr(,"a1_error_metadata")
#> attr(,"a1_error_metadata")$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$status
#> [1] "approximate"
#> 
#> attr(,"a1_error_metadata")$verified
#> [1] FALSE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature
#> attr(,"a1_error_metadata")$exact_moment_quadrature$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$engine
#> [1] "gauss-laguerre"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_verification
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$M_verification_required
#> [1] 160
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_available
#> [1] TRUE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$status
#> [1] "approximate"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$reason
#> [1] "fixed_order_unverified"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_performed
#> [1] FALSE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification_passed
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$absolute_tolerance
#> [1] 1e-10
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$relative_tolerance
#> [1] 1e-08
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$mean_difference
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$variance_difference
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$mean_tolerance
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$variance_tolerance
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$distribution
#> [1] "Gamma(shape, rate)"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$shape
#> [1] 2
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$rate
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$M_verification
#> [1] NA
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$normalization
#> [1] "log-sum-exp"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$normalized_weight_sum
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$schema_version
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$engine
#> [1] "golub-welsch-dense-eigen"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$M_selected
#> [1] 80
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$alpha_param
#> [1] 1
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$cache_key
#> [1] "M=80|alpha=1"
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$cache_hit
#> [1] TRUE
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$zero_normalized_weights
#> [1] 29
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$selected$node_rule$normalized_weight_sum
#> [1] 1
#> 
#> 
#> 
#> attr(,"a1_error_metadata")$exact_moment_quadrature$verification
#> NULL
#> 
#> 
#> attr(,"a1_error_metadata")$proxy
#> attr(,"a1_error_metadata")$proxy$conditional
#> [1] "Poisson(alpha * cJ)"
#> 
#> attr(,"a1_error_metadata")$proxy$marginal
#> [1] "1 + NegativeBinomial(a, b/(b+cJ))"
#> 
#> attr(,"a1_error_metadata")$proxy$cJ
#> [1] 4.47919
#> 
#>