Skip to contents

Generates random samples from the rho = sum w_h^2 distribution by stick-breaking simulation.

Usage

rrho(
  n,
  a,
  b,
  n_sticks = 500L,
  remainder_tol = 1e-10,
  max_sticks = 100000L,
  strict = TRUE
)

Arguments

n

Integer; number of samples to generate.

a

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

b

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

n_sticks

Integer; minimum number of sticks generated for each draw. The sampler continues adaptively when the remaining squared mass is above remainder_tol. Default is 500.

remainder_tol

Numeric; maximum allowed deterministic upper bound on the omitted contribution to rho. Because all ungenerated weights sum to remainder \(R\), their contribution satisfies \(0 <= \sum_{\mathrm{tail}} w_h^2 <= R^2\). Default is 1e-10.

max_sticks

Integer; hard per-draw stick ceiling. Default is 100,000.

strict

Logical; if TRUE (the default), reaching max_sticks before the remainder contract is met raises a typed dpprior_rho_truncation_error. If FALSE, the finite lower approximation is returned with status "approximate" and its deterministic remainder bound in the "rrho_diagnostics" attribute.

Value

Numeric vector of length n. The "rrho_diagnostics" attribute records the adaptive method, status, tolerance, sticks used, and the per-draw deterministic omitted-rho bounds on both ordinary and log scales. When a positive bound is smaller than the floating-point range, its ordinary-scale representation is rounded upward to the smallest positive double rather than silently reported as zero.

Details

Uses the hierarchical representation:

  1. alpha ~ Gamma(a, b)

  2. v_h | alpha ~ Beta(1, alpha) independently

  3. w_1 = v_1, w_h = v_h * prod(1 - v_l) for l < h

  4. Generate at least n_sticks terms and continue until the unallocated remainder \(R\) obeys \(R^2 <= remainder_tol\)

  5. Return the partial sum; its omitted contribution is in \([0,R^2]\)

Thus finite-stick error is controlled draw by draw rather than being silently discarded. Useful for Monte Carlo validation of analytical formulas.

See also

rw1 for w1 random generation

Examples

set.seed(42)
rho_samples <- rrho(1000, a = 2, b = 1)
mean(rho_samples)
#> [1] 0.4205994
mean_rho(a = 2, b = 1)
#> [1] 0.4036526