Control numerical integration of selection likelihoods — set_selection_likelihood

Creates numerical integration settings for the product and best-p-value selection likelihoods used by bselmodel(), bselmodel.mv(), RoBMA(), and RoBMA.mv(). Covariance factors use deterministic sequences of Gauss-Hermite rules. Certified factors of structural rank one through eight use deterministic interval refinement for rank-one product fallbacks and randomized quasi-Monte Carlo for other unresolved factor integrals. General covariance blocks with product weights first try same-diagonal covariance approximations with a common Gaussian factor and independent group factors. Monotone weights permit covariance bounds on the normalizer. Two-sided and other nonmonotone weights require a covariance-perturbation bound; agreement of auxiliary normalizers alone is insufficient. The absolute bound uses the total variation of the step weights. Strictly positive weights also permit a relative Gaussian covariance bound, including reconstructed diagonal and factor rounding, which remains useful for tiny normalizers. Their original full covariance remains authoritative. The covariance error bound, estimated quadrature error, and a bound on Gaussian mass outside the quadrature nodes share one numerical error budget. Interior quadrature error is estimated, so this is not a rigorous probability interval. Unsupported or unresolved cases use randomized quasi-Monte Carlo in the dense conditional-normal representation. Under product weights, a deterministic rung of the factor or cluster quadrature ladder is accepted by a shift-amplitude certificate. Shifting the nodes of one latent axis by a fraction of the node gap and reweighting leaves the integral unchanged, so the change of the rule under half and quarter gap shifts of each axis measures its aliasing error. The rung's successive changes must be within relative_tolerance, the summed shift amplitudes must not exceed half of it, and twice their sum, together with the exterior Gaussian mass, must fit within it. A rung must also resolve the selection edges: the edge width (conditional residual SD divided by the absolute loading) is at least 0.15 times the node gap. When no rung is accepted, rank-one clusters, tensors of rank at most two and depth-one forests try the average of the finest rung over four shifts per axis before sparse grids and the randomized fallback. That average asks the edges for at least one eighth of the node gap, with a separate penalty for each axis's row count, and is used only when relative_tolerance is at least 0.005, the tolerance at which its checks were calibrated. For product weights with selection transitions, the ladder and this refinement tier use the smaller of relative_tolerance and 0.005 for their acceptance checks: a loose request must not relax the calibration of these deterministic rules. After these checks decline, rank-one product integrals refine bounds over the complete Gaussian line, including both tails. Exact row groups and nonnegative excess weights give interval extrema; the returned midpoint must satisfy the requested relative tolerance against both bounds. Floating arithmetic uses an explicit working-precision allowance, rather than a universal guarantee for library probability evaluations. Insufficient work or unresolved arithmetic is rejected without a QMC retry. Other randomized fallbacks retain the requested tolerance and their empirical diagnostics. The best-p-value rule and the other deterministic routes instead require that the node gap across an interior transition not exceed its conditional residual SD divided by the absolute loading; coupled axes use conservative central node spacing. Unresolved candidates refine or use the existing fallback; the numerical settings and likelihood target are unchanged. These checks are heuristics, not a bound on the total integration error, including reciprocal normalizers. Each fit or post-fit request uses a fixed integration design, so repeated likelihood evaluations are deterministic. One-observation blocks are evaluated analytically.

Post-fit dense two-bin product normalizers can use a bounded interpolation of the standard-normal CDF. Its deterministic error is included in the same numerical error budget, separately from Monte Carlo standard error. Other weights, arithmetic outside its supported range, or insufficient error budgets retain the ordinary scalar primitive. An unresolved interpolation attempt retries the ordinary covariance-envelope calculation before QMC. Fitting and random generation retain their ordinary scalar primitive.

The numerical defaults are package operating settings, not statistical estimands or literature-mandated constants. points_per_scramble, max_points_per_scramble, and scrambles set the fixed computational design, while relative_tolerance sets its diagnostic acceptance rule. Increasing a budget changes only the numerical evaluation of the same likelihood target.

For qCMDE/IWMDE post-fit densities and hypotheses in Gaussian selection models, pass this object as density_control$integration_control, for example density_control = list(integration_control = set_selection_likelihood_control(max_points_per_scramble = 32768)). The default NULL retains the fitted integration settings. An explicit control rebuilds a temporary execution plan from the stored covariance metadata without changing the fitted object or its posterior draws. In fitting and qCMDE/IWMDE likelihood evaluation, the maximum point budget controls factor QMC fallback and the rank-one product refinement work ceiling. For rank-one refinement, max_points_per_scramble and scrambles allow at most their product (with max_points_per_scramble rounded up to an even number) times the active exact row-group count times the positive excess-weight bin count in scalar probability evaluations. This includes latent interval masses and bound preparation. Active intervals are limited to the same product and to the existing native node ceiling. points_per_scramble initializes QMC routes only. Post-fit full-event and partial-event calculations, including z-plots and sensitivity diagnostics, also refine dense QMC evaluations up to that budget. Analytic and deterministic quadrature rules remain unchanged.

Usage

set_selection_likelihood_control(
  points_per_scramble = 512L,
  max_points_per_scramble = 8192L,
  scrambles = 8L,
  relative_tolerance = 0.005,
  seed = 1L
)

Arguments

points_per_scramble: nominal number of integrand evaluations per randomized quasi-Monte Carlo scramble for certified factor fallbacks, and the number of shifted Halton points for general covariance blocks. Factor likelihood-normalizer fallbacks under product weights combine the prior Gaussian with a Gaussian at the tilted mean and covariance of the finest quadrature rung and a copy with 1.5 times its standard deviations, where the quadrature provides them (tensors and depth-one forests with at least three finite rules, exterior mass covered and a positive definite covariance). Otherwise, and for the best-p-value rule, they combine the prior Gaussian with proposals optimized from zero and selection-boundary starting points. Every finite optimized candidate is retained, including candidates that locate the same mode. The total count is divided as evenly as possible across the proposals, and importance weights use the actual allocation. Odd likelihood-normalizer factor budgets use the next even total because the fixed design uses two streams. This count is the initial nested-design budget, not an accuracy guarantee. The integration diagnostics determine whether it is adequate for each evaluated likelihood state. Selected-density projection fallbacks instead retain the prior Gaussian and one proposal optimized from zero.
max_points_per_scramble: largest nested-design budget used for a certified factor fallback, including fitting-time refinement, or a post-fit Gaussian event calculation, when the initial points_per_scramble budget fails its integration diagnostic. Refinement doubles the point count up to this value while reusing the same fixed design. Rank-one product fallback instead uses this maximum and scrambles as the deterministic scalar-probability work ceiling described above.
scrambles: number of independent randomized shifts used to estimate integration error for certified factors and general covariance blocks. The stopping multiplier of relative_tolerance (2.36) is that of the default eight scrambles and does not change with this number. It also scales the rank-one product deterministic work ceiling.
relative_tolerance: largest accepted successive relative quadrature change for factor integration, or the largest nested-design relative change and 2.36 times the relative Monte Carlo standard error for a factor fallback and for the conditioned-sampling normalizer. Rank-one product interval refinement instead checks the conservative relative radius of the actual returned midpoint against both bounds. General covariance blocks and Gaussian event integrals use 2.36 times the relative Monte Carlo standard error alone. Dense covariance-envelope quadrature additionally requires the sum of covariance-bound width, twice the successive quadrature change, the relative exterior-mass bound, and any post-fit CDF interpolation error to satisfy this tolerance. Factor quadrature also checks exterior Gaussian mass before accepting rule agreement; an unresolved rule sequence uses the existing QMC fallback. The returned deterministic candidate also requires its shift-amplitude certificate under product weights (summed shift amplitudes within half of this tolerance and twice their sum, with the exterior mass, within the tolerance; the four-shift average that follows a rejected ladder is used only for tolerances of at least 0.005) or, on the other routes, selection-transition resolution; coarser comparison rules retain their actual values without an added resolution requirement. The randomized fallback stops when the larger of its relative change and 2.36 times its relative Monte Carlo standard error (the 97.5% quantile of Student's t distribution with seven degrees of freedom, for eight scrambles) is within this tolerance, and a result that does not pass that test is rejected; the reported standard error stays unscaled, and the criterion does not bound the error of the accepted value. Zplot projections additionally check density error relative to the peak within each posterior draw and raw normalizer-mass diagnostics, including their bounded approximation and omission terms. These diagnostics do not certify the total integrated error of a compressed displayed curve; see as_zplot.brma().
seed: non-negative integer used only to create fixed quasi-Monte Carlo designs. It does not alter R's random-number state.

Value

A RoBMA_selection_likelihood_control object.