Pathfinder inits fail with "Assertion on 'weights' failed: Element <i> is not >= 0" #1323

@jgabry, I'd like your opinion on what cmdstanr should do here

Problem

$sample(init = pathfinder_fit) can fail with:

Assertion on 'weights' failed: Element 163 is not >= 0.

The error comes from posterior::resample_draws() in process_init_approx(). This happened in a cluster run that used Pathfinder inits with num_paths = 10, single_path_draws = 40, draws = 4, psis_resample = FALSE. The Pathfinder seed was fixed, so the failure was deterministic, and a retry gave the same error. Versions: cmdstanr 0.9.0.9000, posterior 1.7.0. The code path is the same in current master (55342f4).

Cause

process_init_approx() (https://github.com/stan-dev/cmdstanr/blob/55342f4f6c/R/args.R#L1381-L1384) smooths the importance weights on the natural scale:

weights <- posterior::pareto_smooth(
  exp(log_weights - max(log_weights)), tail = "right", r_eff=1, return_k=FALSE)

When most Pathfinder paths are poor, lp__ - lp_approx__ is more than about 745 below the maximum for most draws. exp() of these values underflows to exactly 0, but the log weights stay distinct, so the smoothing branch is used. When more than ndraws - ndraws_tail weights are 0, the tail cutoff is 0 and ties with tail values. posterior::ps_tail() then moves the GPD location to 0 - .Machine$double.eps, and some smoothed weights are about -2e-16. resample_draws() rejects them.

Minimal example without Stan:

library(posterior)
w <- c(rep(0, 350), exp(-seq(0, 60, length.out = 50)))
ws <- pareto_smooth(w, tail = "right", r_eff = 1, return_k = FALSE, verbose = FALSE)
sum(ws < 0)  # 7 with posterior 1.7.1
x <- as_draws_df(data.frame(theta = seq_along(ws)))
resample_draws(x, ndraws = 4, weights = ws, method = "simple_no_replace")
#> Error: Assertion on 'weights' failed: Element 341 is not >= 0.

stan-dev/posterior#468 fixes the cause in ps_tail(): it removes the eps offset, so the smoothed tail cannot go below the cutoff. With that fix, the example gives no negative weights, and the natural-scale smoothing gives the same result as the log-scale smoothing (option 2 below).

Options for cmdstanr

cmdstanr requires posterior (>= 1.5.0), so users with an older posterior still get the error after the posterior fix. Possible changes:

1. Clamp after smoothing (one line):

   weights <- pmax(weights, 0)

   This works with every posterior version and adds no warning. The rowsum() collapse that follows then sums non-negative values.

2. Smooth on the log scale, which is the posterior API for importance weights:

   weights <- exp(posterior::pareto_smooth(
     log_weights - max(log_weights), are_log_weights = TRUE, r_eff = 1, return_k = FALSE))

   This does not have the defect in any posterior version, because exp() of the result is never negative. But log_weights can contain -Inf (L1370). For input with -Inf, pareto_smooth() does not smooth and gives the warning "Input contains infinite or NA values ...". To prevent the new warning, cmdstanr would have to smooth only the finite entries.

3. No change in cmdstanr. Require the posterior version that contains the fix, after it is released.

I guess my preference is option 1, because it is the smallest change and it protects users with any posterior version. Which option do you prefer?