prior_truncated() now samples and normalizes correctly when the lower bound is far in the upper tail. The truncated mass and the inverse-CDF uniform came from the lower-tail CDF alone, so for prior_truncated(prior_normal(0, 1), lower = 8) the CDF at the bound rounded to 1 and 5000 draws held 7 distinct values, some below the bound and 440 of them Inf. log_prob() used a wrong normaliser (1.221 instead of 1.289 at 8.1), and prior_truncated(prior_exponential(1), lower = 40) errored with "no probability mass". new_marginal() now works on the log scale of whichever tail holds the bounds (the survival function when the lower bound is past the median), and marginal_sample() draws on that same tail, clamped to the bounds. Left-tail and central truncation give the same results as before up to rounding (#385) (#386).seed = no longer reseeds the caller's global R RNG. npe(), nle(), nre(), sbc(), tarp(), c2st(), simulate_for_sbi() and npe_sequential() each called a bare set.seed(seed), so the user's next runif() depended on seed. A user running npe(..., seed = k) or sbc(..., seed = k) in a loop got downstream draws that repeated or correlated across iterations. A new internal local_seed() seeds the function and restores the caller's .Random.seed (or removes it if there was none) when the function exits, and all eight entry points use it. Seeded results are unchanged. npe_sequential() still derives its per-round seeds from the stream it seeded, inside the restored scope. Torch's generator was already restored by train_restarts() and c2st() (#381) (#383).sample() now returns n draws from a bounded prior at low acceptance. Each rejection-sampling round in sample.nsbi_posterior() asked de_sample() only for the shortfall, so at acceptance a the shortfall shrank by a factor of 1 - a per round and the 100-round cap ran out near 50 * n attempts. With prior_uniform(0, 1) and x_obs = 1.3 (about 1% acceptance), sample(p, 1000) returned 678 rows with a warning, and map_estimate(), sbc(), tarp() and npe_sequential() failed on the short draw. The first round still asks for n; later rounds now size the batch as 1.2 * shortfall / acceptance, with the running acceptance floored at 1% and the batch capped at 100 * n, and call de_sample() in chunks of at most 1e5 rows to bound memory. The warning when max_sampling_batches runs out is unchanged (#380) (#382).npe(), nle() and nre() now error when only one of theta and x is supplied. prepare_simulations() uses precomputed simulations only when both are given, so a lone theta or x was dropped without a word while the simulator ran its full n_simulations budget. The documented contract is that simulator and n_simulations are ignored only when the pair is supplied. A new internal check_theta_x_pair() runs in the early argument validation of all three fitters, so the error ("theta and x must be supplied together") fires before any simulation (#377) (#379).sbc()'s uniformity p-value no longer rejects exactly uniform ranks when the rank bins are unequal. The L + 1 possible rank values 0..L were cut into 20 equal-width intervals, which hold 2 or 3 values each at L = 50 and 51 or 50 at L = 1000, but chisq.test() was called with its default equal cell probabilities. Exactly uniform ranks were rejected 100% of the time at L = 50, 10% at L = 100 and 7% at L = 1000, so a calibrated fit got a small p-value at small n_posterior_samples or large n_sbc. A new internal sbc_rank_bins() splits the rank values into near-equal groups and returns each bin's share, which sbc() now passes as p. plot_sbc() uses the same bins, so its expected-count line and band follow each bin's share instead of the constant n_sbc / bins (#376) (#378).linear_gaussian now warns when there are too few simulations for its regression. fit_linear_gaussian() regresses theta on x with an intercept, so when the number of training rows is at most dim_x + 1 the fit interpolates the data, the residuals are exactly zero, and the posterior covariance reduces to the ridge. The result was a near-delta posterior with huge log-densities and no signal to the user. The fit now warns, naming the row and coefficient counts and the two fixes (simulate more, or raise ridge) (#374) (#375).nsbi_batch_apply() now cancels still-pending futures when a batch's simulator call errors, instead of abandoning them. The parallel collection loop in R/parallel.R calls future::value() on each newly-resolved future; when a batch's simulator call had errored, future::value() re-throws and the error propagated straight out of nsbi_batch_apply()/run_simulator(), leaving any other already-dispatched, not-yet-collected futures running on their workers with their results discarded. This is the dispatch loop behind every simulator call in the package under a future plan (npe(), simulate_for_sbi(), npe_sequential(), sbc(), tarp(), posterior_predictive()), so a diverging simulator or a bad parameter draw wasted worker compute and could leave multisession workers busy after the caller had already handled the error. The collection loop now calls future::cancel() (best-effort, wrapped in its own tryCatch since cancellation is backend-dependent and can itself fail) on every future still pending before re-throwing the original error; the error itself still propagates unchanged (#372) (#373).train_conditional_de() no longer crashes when an epoch's validation loss comes back NaN. Every epoch called scheduler$step(val_loss) unconditionally, but lr_reduce_on_plateau()'s internal .is_better() comparison has no NaN guard, and R's NaN < x evaluates to NA rather than Python's False, so a NaN val_loss -- an exploding-gradient restart, say -- threw "missing value where TRUE/FALSE needed" from inside the scheduler and aborted the whole call, losing every earlier restart's progress. A NaN epoch was already treated as "no improvement" for best-value tracking a few lines above; the scheduler step now gets the same is.finite(val_loss) guard, so it just skips that epoch's LR decay instead of crashing (#370) (#371).npe_sequential() no longer silently biases a round's truncation threshold when the round-2-or-later reference draw comes back short. The r > 1 branch of the round loop calls sample.nsbi_posterior() to draw n_truncation_samples reference samples and takes their epsilon-quantile log-density as that round's acceptance threshold, but sample.nsbi_posterior() only warns, and returns fewer rows than asked for, when a bounded prior's density estimator leaks mass outside the prior support faster than rejection sampling can keep up (capped by max_sampling_batches). A short reference sample fed straight into stats::quantile() computed a threshold from fewer draws than requested, with no round-specific signal to say so -- the same failure mode diagnostic_draws() already guards against in sbc()/tarp(). npe_sequential() now stop()s with a round-specific message naming the round, the shortfall, and the fix (train on more simulations, or lower n_truncation_samples), instead of proceeding on a biased threshold (#368) (#369).tarp()'s ECP curve no longer biases coverage down at points where a trial's coverage value lands exactly on a level's grid boundary. Each trial's coverage value f is count / n_posterior_samples for an integer count, so f lands exactly on one of the 0.05-spaced ECP levels whenever n_posterior_samples is a multiple of 20 -- and, at the top of the range, f == 1 exactly whenever every posterior draw lands closer to the reference point than the truth does, which is not a rare edge case. The curve was computed with a strict mean(f < a), which excluded a trial tied exactly on the boundary even though it genuinely belongs to the credible region at that level; expected_coverage() had the identical bug, fixed in #308 with a closed u >= lo & u <= hi. tarp() now matches it with mean(f <= a) (#366) (#367).check_function() now verifies arity, so a wrong-arity density_estimator/classifier fails before the simulation budget is spent. npe(), nle(), nre(), and npe_sequential() all call a user-supplied density_estimator/classifier as f(theta_z, x_z) -- two positional arguments -- but check_function() only checked that the function took at least one, documented as the check for simulator/sample_fn/log_prob_fn, which really do take one. A function with one formal and no ... passed silently and only failed once training called it, with an unnamed "unused argument" error and the whole simulation budget already spent. check_function() takes a new n_args argument (default 1, so simulator/sample_fn/log_prob_fn are unaffected); the four call sites that validate density_estimator/classifier now pass n_args = 2, so a wrong-arity function is caught before the simulator runs and the error names the argument and the expected arity (#361) (#363).plot_posterior_predictive() now names x_obs and its expected length on a mismatch, instead of a bare, unnamed stopifnot() error. Every other shape/length check in the package goes through a helper that says what argument was wrong and why; this was the one shape check still using a raw stopifnot(length(x_obs) == d), which errored with "Error: length(x_obs) == d is not TRUE" and no indication of which function raised it, what length was expected, or what was received (#362) (#364).sbi head-to-head benchmark (inst/benchmarks/01-04) has been run for the first time, against sbi 0.27.0. gaussian_linear (10 dimensions, MDN and MAF, 10k simulations) passes the roadmap's C2ST <= 0.60 bar outright, against both sbi and the exact posterior: MDN 0.57-0.58, MAF 0.53. two_moons (2 dimensions, no closed-form reference) misses it: MDN 0.62-0.66 at every one of 5 observations, MAF on the line at 0.59-0.66, passing 3 of 5. Every posterior mean/sd difference stays small on two_moons regardless (under 0.04 mean, 0.02 sd, against a [-1, 1] prior) -- the gap there is c2st()'s MLP classifier catching a difference in the shape of two independently fit posteriors of a bimodal target, with no reference to check either against, not a large or obviously wrong fit on either side. Results, the full reading, and the raw comparison CSVs: docs/benchmarks/npe-vs-sbi.md. docs/verification-roadmap.md's status table and "Next actions" are updated to match (#358) (#360).task_sir()'s simulator no longer overshoots to a physically-impossible trajectory for plausible draws from its own prior. The Euler integrator took a single dt = 1 day step per day, with no clamp on S/I inside the loop: for beta draws in the tail of prior_lognormal(meanlog = log(0.4), sdlog = 0.5), newinf = beta * S * I / N could exceed the current S, driving S negative and I above the population N, then oscillating. This happened for roughly 0.7-1% of draws from the task's own prior at the default N = 1e6, days = 160 -- not freak outliers. The bug was silent downstream: the simulator's own output is clamped to [0, 1] before stats::rbinom(), so a contaminated draw produced a finite, plausible-looking observation indistinguishable from a legitimate one, and drop_failed_sims() only screens for NA/NaN/Inf. The integrator now takes 20 Euler sub-steps per day (sir_trajectory(), a new internal helper in R/tasks.R), which shrinks each step's increment enough that newinf/newrec stay small relative to S/I across the prior's whole support; each sub-step additionally clamps newinf/newrec to not exceed the current S/I, as a second line of defense for parameter draws extreme enough to overshoot even a sub-stepped step. task_sir()'s simulator is unchanged in its public contract -- same prior, same N/days/n_points/n_obs_draws arguments, same output shape -- and sir_trajectory() also makes the full daily S/I/R path inspectable, which the strengthened regression test in test-tasks.R uses to draw 2000 theta from the prior and check S/I stay within [0, N] throughout, rather than only checking the already-clamped [0, 1] output (#357) (#359).nre() no longer crashes with a cryptic sprintf error when n_simulations is fractional and the classifier is a neural one. nre() builds n_hint, a row-count hint used to validate the training controls before any simulation runs, straight from n_simulations, but only checked that it was finite and at least 1, not that it was a whole number. A fractional value like n_simulations = 15.5 flowed into check_train_controls(n = n_hint, ...), whose n-dependent stop() branches format n with sprintf's "%d" -- which errors outright on a non-integer double ("invalid format '%d'; use format %f, %e, %g or %a for numeric objects") instead of naming the real problem. classifier = "logistic" never hit this: its closed-form fit takes a looser min_val_rows floor that n_simulations = 15.5 clears, so only the neural classifiers ("resnet", the default, plus "mlp" and "linear") reached the crash. n_hint now mirrors check_count()'s own whole-number test, so a fractional n_simulations falls through to n_hint = NULL for every classifier and defers to prepare_simulations()'s existing check_count() error, the same friendly message "logistic" already gave (#355) (#356).sample() and log_prob() on an nsbi_posterior now reject a zero-row obs/x instead of crashing with a raw "subscript out of bounds". check_matrix() validated column count but enforced no minimum row count, so a zero-row observation (real data filtered down to nothing, then passed straight through as x_obs) sailed past it and reached resolve_obs()'s x[1, , drop = FALSE], which R fails on with an unnamed indexing error rather than anything pointing at the argument. check_matrix() now takes an opt-in min_rows argument (default 0, so every other caller -- within_support(), theta checks, and the existing zero-row test for check_matrix() itself -- is unaffected); resolve_obs() passes min_rows = 1, so sample(), log_prob(), and a zero-row x_obs stored on the posterior at construction all now name the argument and say "must have at least 1 row" instead of crashing (#349) (#353).c2st() now accepts n_folds equal to the smaller sample set's size, the valid boundary, instead of rejecting it. The validation at the top of c2st() read if (n_folds >= n_each) stop(...), which rejected n_folds == n_each even though c2st_stratified_folds() handles that case correctly: folds are cut within x and within y separately via rep_len(seq_len(n_folds), n), and n_folds == n_each still gives every fold exactly one row of the smaller class, so no fold's test set is empty and no fold AUC comes back NA. Only n_folds past that point leaves some fold indices with zero rows of a class. The check is now n_folds > n_each, matching what c2st_stratified_folds()'s own docstring already documented as the real constraint (#350) (#354).npe()/nle() no longer silently fit a garbage model on zero-row precomputed theta/x. prepare_simulations() (R/npe.R), shared by npe(), nle() and nre(), checks n_simulations >= 2 on the simulator path but had no equivalent floor on the precomputed-data path: npe(prior, theta = matrix(numeric(0), 0, 2), x = matrix(numeric(0), 0, 2), density_estimator = "linear_gaussian") ran to completion, with fit_standardizer() taking colMeans()/sd() of zero rows and the resulting NaN center surviving standardization unnoticed. prepare_simulations() now runs theta/x through a new check_min_rows() (R/check.R) before fitting either standardizer, naming whichever argument is empty. nre() was already non-silent about the fully-empty case, via its own pre-flight check_train_controls(n = nrow(theta), ...) call, just with a message about validation_fraction rather than the empty argument; it now also gets check_min_rows()'s clearer message for the case that check never covered, a precomputed x empty while theta is not. A single precomputed row is untouched by this fix and still trains, as before, with warn_constant_columns()'s existing warning on the way past standardization -- that warning's message no longer hardcodes "has one row" for what could be zero rows too, and now reports the actual count (#348) (#351).npe_sequential() now rejects a non-logical standardize before round 1 simulates, instead of after burning its whole simulation budget. npe_sequential() pre-checks several of the npe() arguments passed through ... -- n_components/n_transforms/hidden/n_bins/tail_bound, the training controls, device, embedding_net, density_estimator -- for exactly this reason: a typo should fail immediately rather than after sample_prior() and the user's simulator have already run round 1's full budget. standardize is validated inside npe() itself with check_flag(), but that call only runs once npe_sequential() reaches the npe() call at the bottom of round 1's loop, so it was left out of the pre-check battery that verbose and the rest already get. npe_sequential() now validates standardize with check_flag() alongside the other npe_arg(...) pre-checks, the same pattern verbose already uses (#346) (#347).log_prob() on an NLE/NRE MCMC posterior now rejects a non-logical normalize instead of misreading it. log_prob.nsbi_mcmc_posterior() tested normalize with !missing(normalize) && isTRUE(normalize), and isTRUE() only recognizes the literal value TRUE: normalize = 1 or normalize = "yes" silently took the "unnormalized, no warning" branch instead of raising the "normalize must be TRUE/FALSE" error its sibling log_prob.nsbi_posterior() already raises for the same typo. This is the same bug class already fixed for normalize on log_prob.nsbi_posterior() and standardize in npe()/nle()/nre() (#340) (#342), and, before that, for sum_iid (#321) (#326), z_score in c2st() (#327) (#330), model in the Stan export functions (#328) (#331), refresh in mcmc_draws() (#329) (#332), verbose across the package (#336) (#337), and the save_npe() estimator-class check (#334) (#335). log_prob.nsbi_mcmc_posterior() now validates any explicitly supplied normalize with the existing check_flag() helper (R/check.R) before deciding whether to warn, so a typo is caught at the call that made it and named in the error; the default (normalize left unset, which stays silently unnormalized) is unaffected, since mcmc_log_prob()'s "no normalizing constant" warning was only ever meant to fire when the caller explicitly asks for normalization (#344) (#345).sample() and log_prob() on a posterior now re-check a torch-backed fit's network on every call, not only at posterior() construction. check_fit_alive() exists to fail at the door rather than three calls later, but posterior.nsbi_npe() and mcmc_posterior() (the NLE/NRE path) only ran it once, when the posterior object was built. A fit's torch network can go stale on the same object afterward -- typically from saveRDS()/readRDS() across R sessions -- and every downstream call reached the dead pointer directly, surfacing torch's raw "external pointer is not valid" instead of this package's message pointing at save_npe()/load_npe(). sample.nsbi_posterior(), log_prob.nsbi_posterior(), sample.nsbi_mcmc_posterior() and log_prob.nsbi_mcmc_posterior() all resolve their observation through the shared resolve_obs() helper, so the check now runs there instead of being pasted at each call site; map_estimate() gets the same protection for free since it calls sample()/log_prob() internally (#341) (#343).log_prob() on an nsbi_posterior and npe()/nle()/nre() now reject a non-logical normalize/standardize instead of misreading it. log_prob.nsbi_posterior() branched on normalize through if (normalize && bounded): && coerces a numeric right-hand side silently, so normalize = 1/0 worked by luck, while normalize = NA raised base R's unnamed "missing value where TRUE/FALSE needed" and normalize = "yes" raised "invalid 'x' type in 'x && y'" -- neither pointed at the argument. prepare_simulations()'s if (standardize) had the same problem, and it runs inside npe()/nle()/nre() only after the simulator has already spent the simulation budget. This is the same bug class already fixed for verbose (#336) (#337), sum_iid (#321) (#326), z_score in c2st() (#327) (#330), model in the Stan export functions (#328) (#331), refresh in mcmc_draws() (#329) (#332), and the save_npe() estimator-class check (#334) (#335). Both arguments now validate with the existing check_flag() helper (R/check.R) -- normalize at the top of log_prob.nsbi_posterior(), standardize in npe()/nle()/nre() before the simulator runs -- so a typo is caught at the call that made it and named in the error (#340) (#342).stan_code() and stan_data() now give the right error for an improper uniform prior nested inside prior_independent(), instead of claiming the whole prior is arbitrary R code. check_finite_uniform_bounds() correctly detects an infinite-bound prior_uniform() and points the caller at prior_truncated(), but stan_data() and stan_prior_blocks() only called it inside their identical(prior$type, "uniform") branch. When the improper uniform is a component of prior_independent(), the composed prior cannot use its marginals fast path (an improper marginal has no CDF to restate), so it falls back to closures with prior$type == "independent" and prior$params$marginals == NULL. stan_prior_blocks() then fell into its generic is.null(marginals) branch and reported "arbitrary R code with no Stan counterpart" -- a message written for prior_custom(), not for a named family with an unbounded side. is_improper_uniform_prior() (R/prior.R) already recurses into prior_independent()'s components for exactly this case; nothing called it before the composed-prior path gave up. Both functions now call check_finite_uniform_bounds() unconditionally, before branching on prior$type, so the recursive check runs for every prior shape and the composed case gets the same "improper distribution... prior_truncated()" message a bare prior_uniform(-Inf, Inf) already got. A properly bounded prior_independent() is unaffected: is_improper_uniform_prior() returns FALSE for it, same as it always has for "normal" and the other named families (#338) (#339).npe(), nle(), nre(), npe_sequential(), simulate_for_sbi(), sample() on an NLE/NRE MCMC posterior, and train_conditional_de() now reject a non-logical verbose instead of silently staying quiet. verbose_cat() branched on verbose through if (isTRUE(verbose)) cat(...), and isTRUE() only recognizes the literal value TRUE: verbose = 1 or verbose = "yes" took the "stay quiet" branch with no error, so a caller trying to turn on progress output got silence instead, with nothing said about the typo. This is the same bug class already fixed for sum_iid (#321) (#326), z_score in c2st() (#327) (#330), model in the Stan export functions (#328) (#331), refresh in mcmc_draws() (#329) (#332), and the save_npe() estimator-class check (#334) (#335). Every public entry point that takes a verbose argument now validates it with the existing check_flag() helper (R/check.R) before it is used or forwarded, so a typo is caught at the call that made it rather than read off unexplained silence (#336) (#337).save_npe() now rejects a torch-backed fit whose estimator class it cannot rebuild, instead of writing a file load_npe() can never read back. npe()/nle()/nre() accept a custom density_estimator/classifier function as a documented extension point, and that function can return any S3-classed object wrapping a torch module. save_npe() only checked is.null(fit$de$net) before serializing weights, so it wrote a bundle for a custom torch-backed estimator with no complaint; load_npe() then called de_rebuild_net(), whose switch() on class(de)[1L] only recognizes the four built-in estimator classes and errored for anything else, so the round trip failed only at load time, potentially in a later session after the original fit was gone. save_npe() now checks the estimator's class against the same list de_rebuild_net() switches on (de_rebuildable_classes(), the one place both now read) before doing any work, and errors immediately naming the offending class and what is supported, mirroring stan_pack()'s existing check for the Stan export path. A "linear_gaussian" fit is unaffected: it carries no torch net, so it never reaches this check (#334) (#335).sample() on an NLE posterior now rejects a non-logical refresh instead of silently serving stale cached draws. mcmc_draws() branched on refresh through if (!isTRUE(refresh) && ...), and isTRUE() only recognizes the literal value TRUE: refresh = 1 or refresh = "yes" took the "serve the cache" branch with no error, so a caller trying to force a new MCMC run got the old cached draws back instead, with nothing said about it. This is the same bug class already fixed for sum_iid (#321) (#326), z_score in c2st() (#327) (#330), and model in the Stan export functions (#328) (#331). mcmc_draws() now validates refresh with the existing check_flag() helper (R/check.R), right before it branches on it, so a typo is caught at the call rather than read off stale draws (#329) (#332).stan_code(), write_stan_model() and stan_data() now reject a non-logical model instead of silently mishandling the Stan export. All three branched on model through an unguarded isTRUE(model) (or !isTRUE(model)), and isTRUE() only recognizes the literal value TRUE: model = "TRUE" or model = 1 took the "functions block only" branch in stan_code()/write_stan_model(), or the "x_obs is optional" branch in stan_data(), with no error, instead of the runnable model and the required-x_obs check the model = TRUE default promises. This is the same bug class already fixed for sum_iid (#321) (#326) and for z_score in c2st() (#327) (#330). All three functions now validate model with the existing check_flag() helper (R/check.R), right before they branch on it, so a typo is caught at the call rather than read off a silently wrong export (#328) (#331).c2st() now rejects a non-logical z_score instead of silently disabling standardization. c2st() branched on z_score through if (isTRUE(z_score)), and isTRUE() only recognizes the literal value TRUE: z_score = 1 or z_score = "TRUE" took the "don't z-score" branch with no error, training the classifier on the raw scale of x and y instead of the sbibm-matching standardized scale the default promises. c2st() now validates z_score with the check_flag() helper added for the same bug class in sum_iid (#321) (#326), so a typo is caught at the call rather than showing up as an unexplained accuracy number (#327) (#330).log_lik() and log_ratio() now reject a non-logical sum_iid instead of silently changing what they return. Both functions branch on sum_iid through surrogate_score()'s if (!isTRUE(sum_iid)), and isTRUE() only recognizes the literal value TRUE: sum_iid = "yes" or sum_iid = 1 took the "don't sum" branch with no error, returning an n_theta x n_obs matrix instead of the documented per-theta vector. surrogate_score() now validates sum_iid with a new check_flag() helper (R/check.R), matching how max_batch is already checked in the same function, so a typo is caught at the call rather than read off a wrong-shaped result (#321) (#326).width argument now rejects a length that doesn't match the number of parameters, instead of being silently recycled or truncated. slice_sample_run() fed width straight into rep_len(as.numeric(width), dim) under suppressWarnings(), which accepts any length: an over-long vector was truncated and a short one recycled, both with no warning. posterior.nsbi_nle()/posterior.nsbi_nre() document width as one value per parameter dimension (or a single value to broadcast), so a mis-sized vector silently sized the initial slice interval for the wrong coordinate, the only symptom being a slow- or badly-mixing chain with nothing pointing back to the cause. This is the same class of bug #288/#289/#290/#292 already fixed for theta, x_obs, obs, and within_support()'s bounds, all routed through check_matrix()/check_bound() for exactly this reason; width was the one per-parameter vector in the MCMC path left checking only its values, never its length. check_slice_width_length() now checks length(width) %in% c(1, dim) before rep_len() runs, and errors naming the expected length otherwise (#320) (#325).sample() on an NLE/NRE posterior with seed = ... no longer permanently reseeds the caller's global RNG stream. mcmc_draws() called set.seed(ctl$seed) directly with no save/restore, unlike every other RNG touch point in the package (with_fixed_seed(), set_torch_seed(), rng_streams(), with_rng_stream()). A seeded sample() call therefore left .Random.seed at whatever set.seed(ctl$seed) produced, so any later random draw in the same session became a deterministic function of ctl$seed alone regardless of the caller's own RNG state -- the same bug class #272/#274 fixed for surrogate_potential()'s prior probe and #282/#283 fixed for log_prob()'s acceptance-constant draw. map_estimate() and posterior_predictive() call sample() internally and inherited the same leak. The sampler call now runs under with_fixed_seed(ctl$seed, ...), which parks R's RNG at the seed only for that call and restores the caller's prior state afterward, so seed still makes the run reproducible without leaving a trace on anything the caller does next (#319) (#324).summary() on posterior draws no longer crashes when probs has a single element. summary.nsbi_samples() built its quantile columns with t(apply(m, 2, stats::quantile, probs = probs)): apply() returns a length(probs) x ncol(m) matrix in general, which t() fixes up to ncol(m) x length(probs), but simplifies its result to a plain length-ncol(m) vector when length(probs) == 1, so t() produced a 1 x ncol(m) matrix and the following colnames<- (length 1) failed with "length of 'dimnames' [2] not equal to array extent". summary(draws, probs = 0.5) and summary(post, probs = 0.5) therefore crashed on any fit with more than one parameter -- an obvious way to ask for just the median. The dropped dimension is now restored before transposing, the same guard expected_coverage() already used for the analogous vapply() drop, and probs is validated with check_probs() so an out-of-range value errors clearly instead of failing inside quantile() (#318) (#323).npe_sequential(seed = ...) no longer spends part of every round after the first re-simulating round 1's parameter draws. The round loop passed the caller's seed straight through to the per-round npe() call, and npe() calls set.seed(seed) at the top of its own call. For an estimator that consumes no R-level randomness while fitting, which covers linear_gaussian and any caller-supplied density_estimator function, R's RNG was therefore left in the same state at the end of every round: round r + 1 opened from the state round r opened from, sample_prior() handed it the same candidate matrix, and only the acceptance threshold differed. Three rounds of 200 simulations produced 600 training rows of which 516 were unique; four rounds produced 800 rows and 634 unique. Nothing reported this, so a run that asked for reproducibility quietly threw away a large part of its simulation budget and over-weighted whatever region cleared the threshold in two consecutive rounds, the opposite of what the truncation is for. Each round now derives its own seed from the stream npe_sequential() seeds at the top of the call, which keeps the whole run reproducible from the top-level seed and keeps torch's generator seeded for the neural estimators as #215 requires (#317) (#322).stan_code()'s generated MAF likelihood now builds the alpha-clamp matrix-vector product once per transform, instead of once per element. stan_fn_maf() emitted the log-scale head as for (i in 1:P) al_k[i] = fmin(fmax((Walpha_k * prev + balpha_k)[i], -8.0), 8.0);, which recomputes the full Walpha_k * prev + balpha_k product -- O(P x hidden) work -- on every one of its P loop iterations, even though mu_k two lines above builds the analogous product once as a plain vector[P] mu_k = W * prev + b;. That waste compounds across n_transforms stacked MADE transforms (5 by default) and runs on every leapfrog evaluation of every NUTS iteration for any NLE-MAF fit exported to Stan. The product is now built once as vector[P] al_raw_k, and the loop only applies fmin(fmax(al_raw_k[i], -8.0), 8.0); the clamped values are unchanged, this is a performance fix (#314) (#316).stan_code()'s generated MDN likelihood now builds each mixture component once per _sum_lpdf call, instead of once per observation. stan_fn_mdn()'s <name>_from_head() took the MLP's flat output and rebuilt every component's mean and full Cholesky factor -- softplus on the diagonal, a rep_matrix() plus per-entry assignment for each of the K components -- on every call, even inside the i.i.d. loop where head (and so every component) is identical across all N rows: stan_sum_lines() already hoisted the MLP forward pass out of that loop, but the Cholesky assembly downstream of it still ran N times, on every leapfrog step of every NUTS iteration. <name>_head() is now followed by <name>_mu()/<name>_L(), which build the K means and Cholesky factors as Stan arrays, and <name>_from_components() evaluates the mixture log density from those prebuilt arrays; <name>_sum_lpdf() calls _mu()/_L() once in its precompute block (the same mechanism stan_fn_lingauss() already used) and reuses the result across all N rows, and <name>_lpdf() calls the same three functions for its one observation. The generated log density is unchanged (#313) (#315).tarp() now rejects n_tarp below 2 instead of silently returning a degenerate coverage curve. Distances are computed after z-scoring the true parameter draws by their own spread (fit_standardizer(theta_true)), but sd() of a single value is NA, and fit_standardizer()'s guard against zero-spread columns quietly resets that NA to a scale of 1 with no warning, since it was never told which argument theta_true came from. With n_tarp = 1 that leaves theta_z at exactly zero for the one trial. Under the default references = "uniform", the reference box is range(theta_z) per parameter, which collapses to that same zero point, forcing d_truth = 0 and every trial to read as uncovered until the last level -- an ECP curve that jumps from 0 to 1 at the first step, indistinguishable from a badly miscalibrated posterior regardless of fit quality. references = "prior" sidesteps the forced-zero distance but still standardizes by the same degenerate, unscaled std, so its distances are not comparable across parameters either. n_tarp now goes through check_count(n_tarp, "n_tarp", min = 2), the same guard used elsewhere in the package, so both reference modes get a real spread to standardize by (#311) (#312).expected_coverage() no longer excludes SBC ranks that land exactly on the credible-interval boundary. Ranks are integers in {0, ..., L} (L = n_posterior_samples); u <- rank / L approximates the posterior CDF at truth, and a trial was counted as covered only when lo < u < hi, strict on both ends. Whenever lo * L or hi * L is itself an integer -- the default L = 1000 and alpha = 0.9 give lo = 0.05, hi = 0.95, so lo * L = 50 and hi * L = 950 -- a rank landing exactly on that boundary is genuinely inside the central interval but was thrown out, biasing empirical coverage down by about O(1/L) for a perfectly calibrated posterior (899 covered trials read as 0.898 instead of 0.9 out of 1000, say). The bias was systematic, not sampling noise, so it didn't shrink with more SBC trials. The comparison is now closed on both ends (lo <= u <= hi). plot_coverage() calls expected_coverage() directly and picks up the fix with no change of its own (#308) (#310).nre() now rejects batch_size below what its atomic loss needs, instead of training on it silently. minibatches() only ever merges the trailing short batch into the one before it, so a batch_size below the 2-row floor nre_atomic_log_prob() needs is not a one-off: every interior minibatch of every epoch stays that small, and the atomic loss returns a constant zero gradient below 2 rows. With batch_size = 1, that is every minibatch but the one merged trailing batch, so training ran almost entirely on the signal from that one batch, with no error or warning to say so. check_train_controls() now checks batch_size against min_val_rows, the same floor it already enforces on both sides of the train/validation split (#188, #239), and nre() reports it before simulating rather than after (#307) (#309).npe_sequential() now validates embedding_net before round 1 spends its simulation budget. embedding_net reaches round 1 through ..., alongside n_bins, device, and the rest of npe()'s estimator/training-control arguments, and #251/#262 already moved those checks up front so a bad value fails before prepare_simulations() spends round 1's draws. embedding_net was left out of that pass: npe_sequential(prior, simulator, n_rounds = 3, n_simulations = 500, embedding_net = list(bogus = TRUE)) ran the full round-1 simulation budget and only then failed inside the round-1 npe() call with `embedding_net` must be built with embedding_mlp(). The same check npe() already runs now also runs up front in npe_sequential(), so the error is identical but arrives before any simulation happens (#301) (#305).npe()/nre() now warn when embedding_net is supplied alongside a function-valued density_estimator/classifier. The existing check only fired for identical(density_estimator, "linear_gaussian"), which is FALSE for a function value, so a caller-supplied fitter passed to npe(..., density_estimator = my_fitter, embedding_net = embedding_mlp(4)) silently dropped the embedding: fit_density_estimator()/fit_ratio_estimator() forward only theta/x to a custom function, never embedding_net. Both entry points now warn in that case too, symmetric with the "linear_gaussian"/"logistic" cases (#300) (#304).embedding_mlp() no longer silently truncates a non-integer output_dim or hidden width. Its manual checks tested output_dim against < 1 but never against trunc(), so a fractional value passed the guard and was silently floored by as.integer() two lines below; hidden was coerced with as.integer() before its own check ran, so a fractional entry was truncated first and the check never saw the original number. embedding_mlp(output_dim = 2.7, hidden = c(10.9, 5.2)) used to return output_dim = 2L, hidden = c(10L, 5L) with no warning. Both arguments now go through check_count()/check_counts(), the same helpers every other size argument in the package uses, so a fractional value errors instead of being rounded away. An empty hidden still means a single linear map to output_dim, unchanged (#302) (#303).stan_data() now takes a model argument and requires x_obs when it is TRUE. stan_code()'s default (model = TRUE) always emits a data block that declares N and x as required, but stan_data() made x_obs optional and, when it was omitted, returned a list with neither field. Pairing that program with that data list compiled fine in R and then failed deep inside cmdstanr/rstan with an opaque "variable does not exist" error instead of a message from this package. stan_data() now mirrors stan_code()'s model argument (default TRUE) and errors up front, naming x_obs, when model = TRUE and x_obs is NULL; pass model = FALSE to build a data list for a functions-only export, which has no N/x to fill (#298) (#299).npe()/nle()/nre() no longer guess the layout of a flattened multi-parameter theta on the pre-computed theta =/x = path. prepare_simulations() passed a bare theta vector straight to as_theta_matrix(), which reshapes anything that isn't a single row with matrix(theta, ncol = d, byrow = TRUE) -- silently assuming row-major order. When length(theta) happened to be an exact multiple of prior$dim, that assumption was never checked against anything: theta_flat <- as.vector(theta_matrix), the ordinary column-major way to flatten an n x d matrix in R, came back scrambled, and training proceeded on the wrong parameter values with no warning or error. A single simulation (length(theta) == prior$dim) stays unambiguous and keeps working as before; for a multi-parameter prior, any other bare vector now errors and asks for a matrix or data frame instead of guessing (#291) (#295).
within_support() no longer silently reshapes a wrong-length theta into several parameter sets. It handed theta to as_theta_matrix(), which reshapes whatever length it is given rather than checking it matches prior$dim -- so within_support(prior_uniform(low = c(-2, -2), high = c(2, 2)), c(0.1, 0.2, 0.3)) recycled the length-3 vector into a 2-row matrix with a cryptic recycling warning, and within_support(prior, c(0.1, 0.2, 5.0, 0.3)) (length 4, a clean multiple of dim = 2) silently became 2 rows with no warning at all. This is the same bug class #289 fixed at five other call sites; within_support() was scoped out of that PR for separate consideration. It now goes through check_matrix() like the rest: a bare vector is read as a single row and must have exactly prior$dim entries, and a mismatched length errors instead of recycling (#292) (#293).
posterior() on an npe fit no longer silently reshapes a wrong-length x_obs. posterior.nsbi_npe() handed x_obs to as_theta_matrix(), which reshapes whatever length it is given rather than checking it matches fit$dim_x -- the exact bug class #288/#289 fixed at five other call sites, in a spot that fix missed. posterior(fit, x_obs = c(10, 20, 30, 40)) against a dim_x = 2 fit silently became a 2-row x_obs with no warning at all, since length 4 is a clean multiple of 2. mcmc_posterior(), the sibling path for NLE/NRE, already used check_matrix() for this; posterior.nsbi_npe() now does too, matching the error posterior() on an NLE/NRE fit already gives for the same shape mistake (#290) (#294).
A wrong-length theta/x_obs no longer gets silently reshaped into several parameter sets or observations. log_prob.nsbi_posterior()'s theta, resolve_obs()'s x/obs (backing sample(), log_prob() and map_estimate()), mcmc_posterior()'s x_obs (shared by posterior.nsbi_nle()/posterior.nsbi_nre()), likelihood_fn()'s x_obs, and stan_data()'s x_obs all handed their checked value to as_theta_matrix(), which reshapes whatever length it is given rather than checking it matches the fit's dimension -- so log_prob(post, theta = c(0.1, 0.2, 0.3, 0.4)) against a two-parameter fit silently became a 2-row theta and returned a plausible-looking length-2 answer instead of an error. This was worse for NLE/NRE posteriors: resolve_x_iid() treats every row as an independent observation with no row-count check to catch the reshape at all, so sample(nle_post, obs = c(0.1, 0.2, 0.9, 0.9)) on a dim_x = 2 fit silently conditioned on 2 fabricated observations. All five call sites now go through check_matrix(), which check_x_obs() already used for npe_sequential(): a bare vector is read as a single row and must have exactly the fit's dimension, and a mismatched length errors instead of recycling (#288) (#289).
npe_sequential() no longer reports round 1's acceptance as a hardcoded 1.00 when the simulator fails on some prior draws. The main loop initialized acceptance <- 1 for round 1 and only recomputed it from what drop_failed_sims() actually kept when r > 1L, even though drop_failed_sims() runs -- and can shrink theta_new/x_new by dropping non-finite rows -- for round 1 too. A simulator that routinely fails on some parameter draws (an ODE solver that diverges for certain values, say) then printed a self-contradictory verbose log, e.g. "Round 1/3: 850 new simulations (850 total), proposal acceptance 1.00" after 1000 were requested, and left the wrong value in the returned fit$rounds[[1]]$acceptance and in print.nsbi_snpe()'s acceptance-per-round line. acceptance <- nrow(theta_new) / max(tried, 1L) now runs unconditionally after drop_failed_sims(), with tried set to round 1's full simulation budget since round 1 has no rejection-sampling stage to have already tried fewer draws (#285) (#287).
sample()'s shortfall warning no longer blames prior leakage for an unbounded prior. When a density estimator's draws come up short of the requested count, sample.nsbi_posterior() warned "Only X/Y samples inside prior support ... The estimator is leaking mass outside the prior" regardless of whether the prior was bounded. For a bounded prior (prior_uniform(), or prior_custom() with lower/upper) that wording is correct: the shortfall comes from rejection-sampling draws outside the prior's actual support. For an unbounded prior (prior_normal() and the like) there is no support boundary to leak past -- the shortfall there is the finite-row filter dropping NaN/Inf draws that de_sample() itself produced, and the warning was pointing at the wrong fix. The warning now branches on bounded: the bounded wording is unchanged, and the unbounded wording says the density estimator produced non-finite draws (#284) (#286).log_prob() on a bounded-prior posterior no longer mutates the caller's RNG stream. normalize = TRUE's acceptance-constant estimate draws n_normalization samples from the density estimator via de_sample(), and that draw had no save/restore around it -- the same bug class #274 fixed in surrogate_potential()'s prior probe, just in log_prob.nsbi_posterior() instead. log_prob() reads as a pure evaluation function, unlike sample(), which is documented to consume randomness, so calling it on a bounded prior silently perturbed R's RNG and, for a neural (MAF/NSF/MDN) estimator, torch's global RNG too (de_sample_flow() draws with torch_randn()). The draw now runs under with_fixed_seed(), and, only when the density estimator actually has a torch network (linear_gaussian does not), torch's RNG is saved and restored with set_torch_seed()/torch_set_rng_state(), the same pattern #276 established for a seeded fit or c2st() call (#282) (#283).de_log_prob() on a linear_gaussian estimator no longer crashes on a zero-row theta. Its broadcast check only stretched the single-observation conditional mean mu up to theta's row count when nrow(mu) == 1L && nrow(theta) > 1L; with a zero-row theta that condition is false, so mu stayed a 1-row matrix and dmvnorm_chol(theta, mu, chol) computed theta - mu between a 0-row and a 1-row matrix, failing with "non-conformable arrays". This is the same corner #271 fixed inside dmvnorm_chol() itself (a non-empty mean arriving as a plain vector); here mu is already a matrix, so that fix did not cover it. de_log_prob.nsbi_de_lingauss() now returns numeric(0) for a zero-row theta directly, before lingauss_mean()/dmvnorm_chol() ever run -- the torch estimators (MDN/MAF/NSF) never hit this, since PyTorch's broadcasting already expands a size-1 dimension against size-0 (#279) (#281).c2st(..., classifier = "logistic") no longer errors on a single-column comparison. c2st_logistic_prob() built the training frame as data.frame(y = y_train, x_train) and the prediction frame as data.frame(x_test). data.frame() names a lone unnamed column after the deparsed argument it was given, so an unnamed, one-column x_train/x_test pair landed under different names in the two frames -- x_train in one, x_test in the other. glm() fit a coefficient under whichever name reached the training frame, and predict(fit, newdata = ...) then failed to find it under the prediction frame's name. Both frames now get the same, argument-independent column names before glm()/predict() run, so the fitted formula and newdata line up regardless of column count or the input matrices' own names (#278) (#280).c2st() calls no longer permanently reseed torch's global RNG for the rest of the R session. train_restarts() and c2st()'s MLP path both called torch::torch_manual_seed(seed) directly whenever a caller passed seed, and neither saved nor restored the generator's prior state afterward -- unlike every other RNG touch point in the package, and unlike what #272/#274 already fixed for R's own RNG. Since torch has one global generator with no way to ask for an independent stream, a single seeded npe()/nle()/nre() fit or a seeded c2st() call left every later unseeded torch call in the same session -- another fit, de_sample(), an unseeded c2st() -- drawing from wherever the seeded call left the generator, instead of from a fresh stream. Both call sites now go through set_torch_seed(), which saves torch::torch_get_rng_state() before reseeding so the caller can restore it with on.exit(), the torch analogue of with_fixed_seed() (#275) (#276).log_prob() on an NLE/NRE posterior no longer mutates the caller's RNG stream, and no longer rebuilds surrogate_potential() on every call. mcmc_log_prob() called surrogate_potential() fresh on every log_prob() call, unlike slice_sample_surrogate(), which builds it once per chain and reuses it for every MCMC step; map_estimate()'s optimizer calls log_prob() once per evaluation against the same x_obs, so a Nelder-Mead/BFGS run of a few hundred iterations rebuilt the whole evaluator that many times, re-tensorizing x_obs and never letting an MDN's JIT trace warm up. Inside surrogate_potential(), the prior-log_prob probe (prior$log_prob(sample_prior(prior, 2L))) also drew from the global RNG stream and never restored it, unlike every other RNG touch point in the package -- so set.seed(42); log_prob(post, theta); left .Random.seed different from what set.seed(42) alone would give, breaking reproducibility for any code that called log_prob() in between two seeded random draws. The potential closure is now cached on the posterior (cached_surrogate_potential(), rebuilt only when x_obs/max_batch changes), and the probe runs under with_fixed_seed(), the same save/restore pattern rng_streams() and with_rng_stream() already use (#272) (#274).x/x_obs no longer crashes cross_iid()/mdn_iid_blocks() with a bare base-R error, for every neural estimator and NRE. check_matrix() accepts a 0 x d matrix untouched, and it reaches log_lik(), log_ratio(), and posterior(fit, x_obs = ...) for MAF, NSF, MDN, and NRE alike -- a user who filters an observation set down to nothing, or subsets x_obs by mistake, hit this every time. cross_iid() computed obs_chunk >= 1 and called seq.int(1L, n_obs, by = obs_chunk); with n_obs = 0 that becomes seq.int(1, 0, by = 1), which R rejects with "wrong sign in 'by' argument", naming neither the package nor the actual problem. mdn_iid_blocks() had the identical pattern independently. linear_gaussian's own de_iid_evaluator.nsbi_de_lingauss() already got this right and returns a log-likelihood of 0, the correct empty-product value; cross_iid() and mdn_iid_blocks() now short-circuit on a zero-row theta or x before either seq.int() call runs, since every caller's pre-allocated output is already zero, or already the right shape, and needs no further work. dmvnorm_chol() (R/density_estimator.R) had a related, non-crashing bug in the same corner: recycling a non-empty conditional mean into a zero-row matrix triggered R's "non-empty data for zero-extent matrix" warning on every call, even though the resulting empty matrix was correct; it now builds that shape directly (#271) (#273).c2st() now stratifies its cross-validation folds by class, instead of shuffling one fold assignment over both sample sets together. Its docstring claims to reproduce sbibm/metrics/c2st.py, which reaches sklearn's default StratifiedKFold through cross_val_score() -- but the fold assignment was plain unstratified shuffling (base::sample(rep_len(seq_len(n_folds), n))), so at small sample sizes or with an unlucky shuffle, a fold's test set could land with zero rows of one class. roc_auc() already returns NA_real_ for such a fold (nothing to rank), and mean(aucs) carries no na.rm, so c2st()$auc came back NA with no warning; a 200-seed repro at n_each = 6, n_folds = 5, classifier = "logistic" hit this on 173 seeds. Folds are now assigned within the x-rows and within the y-rows separately (c2st_stratified_folds()), then concatenated, which is what StratifiedKFold guarantees for a two-class target and keeps every fold's test set at least one draw of each class for any n_folds the existing validation allows (#269) (#270).is_improper_uniform_prior() now sees an improper prior_uniform() nested inside prior_independent(), instead of only checking the joint prior's own top-level type. #263 taught surrogate_potential() to detect an infinite-bound prior_uniform() and point at prior_truncated() before probing the prior, but the check only fired when prior$type was literally "uniform". prior_independent(mu = prior_uniform(-Inf, Inf), nu = prior_normal(0, 1)) has type = "independent" regardless of what is inside it, and an improper component forces prior_independent() off its marginals fast path, so the resulting joint prior carried no bounds of its own to inspect either -- is_improper_uniform_prior() returned FALSE, surrogate_potential() fell through to the old probe, and #263's exact misdiagnosis (pointing at prior_custom(..., log_prob_fn = ) instead of prior_truncated()) reappeared one level of composition away. prior_independent() now keeps its component priors themselves in params$components (previously just their type names, discarded once the fast path was skipped), and is_improper_uniform_prior() recurses into them, so the same check covers a prior_uniform() at any depth of prior_independent() nesting without special-casing composite shapes (#267) (#268).posterior()/sample() on an nsbi_nle/nsbi_nre fit now diagnoses an improper prior_uniform() correctly, instead of pointing at the wrong fix. prior_uniform(low = -Inf, high = Inf) (or any infinite bound) is documented and supported specifically so it can be bounded later with prior_truncated(); sample_prior() already refuses to draw from it directly, with a message naming prior_truncated(). surrogate_potential()'s prior-log-density probe called sample_prior() to check the prior works, caught that correct error in a tryCatch(), and turned it into NA_real_ -- so all(is.na(probe)) fired and told the user to rebuild the prior with prior_custom(..., log_prob_fn = ), which is not the problem: prior_uniform()'s log_prob_fn works fine, the prior is just unbounded. R/stan.R already refused this same case for the stan_code()/stan_data() path (#252); surrogate_potential() now runs the same check (factored into is_improper_uniform_prior(), shared by both call sites) before the probe, and reports the actual fix, prior_truncated() (#263) (#266).npe_sequential() now validates density_estimator before round 1 simulates. npe(), nle() and nre() all resolve density_estimator/classifier with match.arg() and call check_torch_for_estimator() before prepare_simulations() runs the simulator (#250) -- but npe_sequential() only ever forwarded density_estimator through ... to the npe() call at the end of round 1, so a typo'd name or a "maf"/"mdn"/"nsf" choice with torch unavailable was only caught after round 1's whole simulation budget was already spent. npe_sequential() now resolves density_estimator and calls check_torch_for_estimator() in the same pre-flight block that already checks n_bins, device, and the rest of round 1's npe() arguments (#251), and passes the resolved value to each round's npe() call rather than re-running match.arg() every round (#262).npe()/nle()/nre() now validate a caller-supplied density_estimator/classifier function before running the simulator, instead of only discovering it is broken at the very end of fit_density_estimator()/fit_ratio_estimator(). Every other argument to these three functions is checked before prepare_simulations() spends the simulation budget (#250), and check_function() already exists for exactly this, used for simulator and for prior_custom()'s sample_fn/log_prob_fn. A caller-supplied estimator or classifier was the one documented extension point that skipped it: npe(prior, simulator, n_simulations = 100000, density_estimator = my_fitter) with a wrong-arity or typo'd my_fitter ran the full 100,000 simulations before failing deep inside training. npe() and nle() now call check_function(density_estimator, "density_estimator", ...), and nre() calls the analogous check on classifier, alongside the other pre-simulation checks; the built-in string options ("maf"/"mdn"/"nsf"/"resnet"/etc.) are untouched (#259) (#261).slice_sample_run() (R/mcmc.R) now coerces a non-finite log-density to -Inf for every candidate it proposes, not just the chains' initial state. The initial state was checked for finiteness up front, but the stepping-out loop and the shrinkage loop each fed log_prob_fn(cand)'s output straight into <=/> comparisons and rowSums()/max.col(), with no guard. log_prob_fn here is surrogate_potential()'s closure, which adds a trained MAF/NSF/MDN/NRE network's output straight into the log-density -- an out-of-distribution theta inside the prior's support but outside the estimator's training distribution can make that NaN, the same failure mode posterior.R and npe_sequential() already guard against (#221/#234/#244). Left unguarded here, a single NaN corrupted an entire batch of candidates at once via NA-tainted rowSums()/max.col(), surfacing as a bare "NAs are not allowed in subscripted assignments" or as NaN silently carried into the retained chain state -- hitting sample() on an nsbi_nle/nsbi_nre posterior, and any sbc()/tarp() run against one, since every trial starts a fresh chain. Both loops now coerce lp[!is.finite(lp)] <- -Inf right after each log_prob_fn() call, the same treatment the package already gives non-finite draws elsewhere (#258) (#260).posterior() on an nsbi_nle/nsbi_nre fit now honors max_batch, instead of silently ignoring it. log_lik() and log_ratio() both chunk the (theta, x) pairs they evaluate according to max_batch, and posterior()'s MCMC path reaches the same evaluator, surrogate_potential() -- but mcmc_posterior() stored a caller's max_batch unchecked into control$dots, and neither slice_sample_surrogate() nor mcmc_log_prob() (which back sample() and log_prob() on the resulting posterior) ever read it, so every MCMC evaluation ran at the default max_batch = 1e5 regardless of what was passed. This bit hardest for nre(), whose ratio classifier has no i.i.d. fast path: every observation costs a forward pass, so there was no way through the public API to shrink a batch and avoid an out-of-memory error on many chains against a large observation set. posterior.nsbi_nle() and posterior.nsbi_nre() now take max_batch as a named argument (default 1e5, matching prior behavior), validated the same way log_lik()/log_ratio() validate it (#230), and threaded through to every surrogate_potential() call the resulting posterior makes (#256) (#257).npe_sequential() now checks round 1's estimator/training-control arguments before round 1 simulates. n_rounds, n_simulations, epsilon, n_truncation_samples and max_proposal_batches were already validated up front, closing this gap for #226, but everything else that round 1 forwards through ... to its internal npe() call -- n_bins, tail_bound, hidden, n_components, n_transforms, batch_size, lr, patience, n_restarts, clip_grad_norm, max_epochs, validation_fraction, device -- was still only checked once that npe() call actually ran, at the end of the round, after prepare_simulations() had already spent round 1's simulation budget. npe_sequential(prior, simulator, x_obs = 0.5, n_rounds = 1, n_simulations = 500, density_estimator = "linear_gaussian", n_bins = 1) ran the simulator 500 times before raising n_bins's "at least 2" error. npe_sequential() now runs npe()'s own check_architecture()/check_train_controls()/check_device_arg() on round 1's ... args before the first prepare_simulations() call, filling in npe()'s own defaults for anything the caller left unset so the two checks cannot drift apart (#251) (#254).
stan_data()/stan_code() now refuse a prior_uniform() with an infinite low/high bound instead of writing it out as nsbi_low = -Inf, nsbi_high = Inf. prior_uniform(low = -Inf, high = Inf) is a legal, if improper, prior -- it exists so prior_truncated() can bound it -- and sample_prior() already refuses to draw from it directly. An nle() fit built from pre-computed theta/x never calls sample_prior(), so the improper prior previously reached stan_data() and stan_prior_blocks() (R/stan.R) untouched: stan_data() returned nsbi_low = -Inf, nsbi_high = Inf, and stan_code() emitted an unconstrained vector[...] theta; with no error from this package, leaving the failure (if any) to surface later inside Stan. Both functions now call a new check_finite_uniform_bounds(), matching the existing prior_custom() rejection in style and message (#252) (#255).
npe()/nle()/nre() now fail before running the simulator when the chosen estimator or device needs torch and torch is unavailable, instead of after. All three validate their cheap arguments -- device, prior, architecture -- before prepare_simulations() runs the simulator, since the simulation budget is the expensive part of a call. require_torch() was only reached deep inside train_restarts() (R/train.R), which runs after prepare_simulations() has already called the simulator n_simulations times, so npe(prior, simulator, n_simulations = 500) (the default density_estimator = "maf") ran the simulator 500 times on a machine without torch before erroring -- the same for nle()'s default estimator and nre()'s default classifier = "resnet". resolve_device()'s CUDA/MPS availability check had the identical ordering problem: device = "cuda" on a non-CUDA machine only failed after simulating. A new check_torch_for_estimator() (R/check.R) runs both checks right after the other cheap-argument checks, for any string density_estimator/classifier that needs torch; density_estimator = "linear_gaussian", classifier = "logistic", and a caller-supplied estimator/classifier function are unaffected (#250) (#253).log_lik()/log_ratio()'s max_batch now bounds memory on the flow/NRE path even when x, not theta, is the large dimension. cross_iid() (R/likelihood.R) chunked theta into blocks sized so theta_chunk * n_obs <= max_batch, but once n_obs alone exceeded max_batch, theta_chunk floored to 1 and that single call still handed score() the entire n_obs-row observation set -- the observation side was never chunked at all, contradicting the "regardless of which side is large" contract restored for the MDN path by #243. This is the path posterior()'s MCMC sampling uses for an NLE or NRE fit, and NLE/NRE's headline use case is conditioning on a large i.i.d. observation set (thousands of trials) without retraining, so this is exactly the shape most likely to spike memory. cross_iid() now blocks theta first as before and chunks observations within each theta block, accumulating into the same flat (theta, x) vector before calling collect() -- no call to score() sees more than max_batch pairs, and every caller's collect(idx, lp) contract is unchanged (#248) (#249).npe_sequential() no longer hands a NaN-filled parameter row to the user's simulator during round 2+'s truncated proposal rejection. Round 2+ candidates are drawn from sample_prior(prior, n_needed) -- the full prior, not the truncated proposal -- which is exactly the out-of-distribution regime where a MAF/NSF density estimator can return NaN from log_prob() (#221). The proposal filter compared that vector to the truncation threshold with a bare >= threshold, so a NaN log-density made keep NA at that row; R's matrix indexing keeps, rather than drops, a row selected by an NA logical index and fills it with NA, so the corrupted row survived into theta_new and reached run_simulator() -- wasting a simulation, or erroring outright for a simulator that checks its input. The loop now guards the comparison with is.finite(lp) & lp >= threshold, the same coercion already used at the log_prob.nsbi_posterior() and mcmc_init_resample()/mcmc_init_proposal() call sites for this exact NA-vs-FALSE problem. fit$rounds[[r]]$acceptance is also now computed after drop_failed_sims() for round 2+, so it reflects proposals the simulator could actually use rather than the pre-simulation count a rejected NaN row used to inflate (#246) (#247).sample() on an unbounded posterior (e.g. prior_normal()) no longer lets a non-finite draw from the density estimator through. sample.nsbi_posterior()'s non-finite-row filter ran only inside if (bounded), since it was added by #234/#236 to work around within_support() returning NA (not FALSE) for a NaN row -- a problem specific to the bounded rejection-sampling path. That left bounded <- !is.null(prior$lower) || !is.null(prior$upper) FALSE for an unbounded prior, the common case in the package's own NPE examples, with no filter at all: a NaN/Inf row that an under-trained MAF/NSF/MDN occasionally produces was rbind'd straight into the returned draws matrix, attr(draws, "acceptance_rate") still reported 1.0, and the corruption propagated silently into summary(), pairplot(), and sbc()/tarp() diagnostics -- or surfaced downstream in map_estimate() as a confusing "theta contains non-finite value" error blaming the seed draw. sample.nsbi_posterior() now drops a non-finite row from de_sample()'s output unconditionally, before the bounded branch's within_support() check runs, so acceptance_rate reflects the drop for every prior (#244) (#245).
log_lik()/log_ratio()'s max_batch now bounds memory for an MDN-based fit even when theta, not x, is the large dimension. For MAF/NSF/linear_gaussian/NRE, cross_iid() (R/likelihood.R) chunks the (theta, x) cross product by theta rows, so max_batch bounds memory regardless of which side is large. The MDN's fast i.i.d. path, mdn_iid_blocks(), only chunked x: it ran mdn_mixture() -- the MLP forward pass and Cholesky assembly -- over the whole of theta in one call before any chunking happened, so scoring a dense theta grid (a profile-likelihood plot, say) against a handful of observations materialized a (n_theta, K, dim_theta, dim_theta) tensor however small max_batch was set. mdn_iid_blocks() now blocks theta first, the same way cross_iid() does, and chunks observations within each block as before; mdn_trace_cache()'s TorchScript shortcut only fires when a single call would cover both dimensions, since a trace recorded at one shape can't stand in for the chunked path (#240) (#243).
nre() no longer silently trains on zero gradient when the training split, not just the validation split, drops to one row. check_train_controls() (R/train.R) enforced min_val_rows (2 for nre()'s atomic contrastive objective, per #188) against n_val only, with no matching floor on n_tr = n - n_val. A large validation_fraction can clear the validation-side floor while leaving n_tr below it -- n_simulations = 4, validation_fraction = 0.75 gives n_val = 3 (passes) and n_tr = 1 (was never checked). train_restarts() then trained on that single row, and nre_atomic_log_prob()'s k < 2L guard -- the same branch #188 fixed for the validation side -- returned a constant zero loss every step: no gradient, no error, training ran to patience epochs and reported a best_val_loss as if it had actually trained. check_train_controls() now also requires n - n_val >= min_val_rows, the same floor already applied to the validation side (#239) (#242).
map_estimate() on a 1-D posterior no longer errors on a prior_custom() with one infinite bound. Its fit$dim_theta == 1L branch chose stats::optim(method = "Brent") whenever prior$lower and prior$upper were both non-NULL, on the assumption that a NULL bound is the only way a prior can be one-sided. prior_custom(dim = 1, lower = -Inf, upper = 5) has both fields set -- one is just non-finite -- so it took the Brent branch anyway, and stats::optim(method = "Brent", lower = -Inf, ...) errored immediately with "'lower' and 'upper' must be finite values" before the search ran at all. The Brent branch now also requires is.finite(prior$lower) && is.finite(prior$upper); a bound that is present but infinite falls through to the existing L-BFGS-B branch, which already handles Inf on the missing side (#238) (#241).n_bins, which could surface as a NaN training loss. rq_spline() (R/nsf.R) turns a softmax over n_bins (K) bins into bin widths via min_bin + (1 - min_bin * K) * softmax, valid only while min_bin * K < 1; check_architecture() enforces n_bins >= 2 but never capped it, so npe(..., density_estimator = "nsf", n_bins = 2000) passed validation with the fixed default min_bin = 1e-3 and 1 - min_bin * K went negative. Once that scale factor is negative, any softmax weight large enough (routine once an autoregressive net starts concentrating mass on one bin) makes its bin width negative, which breaks the rational-quadratic spline's monotonicity and can send its log-determinant term (log() of a non-positive slope ratio) to NaN -- including mid-training, far from n_bins where the actual mistake was made. The identical shape applied to bin heights, built the same way two lines down. rq_spline() now clamps min_bin to min(min_bin, 0.5 / K) before rescaling, so min_bin * K <= 0.5 holds for any n_bins a caller chooses rather than only for n_bins < 1000 (#235) (#237).sample(), log_prob(), and map_estimate() no longer let a NaN draw from the density estimator slip past within_support()'s rejection filter. within_support() (R/prior.R) compares a row against lower/upper with sweep()/rowSums()/==, so a row containing NA/NaN returns NA rather than FALSE. log_prob.nsbi_posterior() already accounts for this on user-supplied theta by running it through check_finite() first (#221), but three call sites test a vector that comes from the density estimator's own output instead, and none of them treated NA as reject: sample.nsbi_posterior()'s rejection filter (draw[within_support(prior, draw), ]) kept an NA-indexed row rather than dropping it -- R's matrix indexing fills it with NA instead -- so a NaN draw from an under-trained MAF/NSF/MDN became an all-NA row counted toward n and returned as a real posterior draw, with no warning; log_prob()'s normalize = TRUE acceptance estimate (mean(within_support(prior, draw))) turned NA and then crashed the following comparison with the unrelated base-R message "missing value where TRUE/FALSE needed"; and map_estimate()'s objective crashed the same way if the optimizer ever proposed a non-finite parameter. All three now coerce within_support()'s NA to FALSE before using it, so a non-finite draw is rejected as leakage the same way an out-of-bounds one is (#234) (#236).log_lik()/log_ratio() on an NLE or NRE fit now reject a non-finite max_batch instead of an opaque rep() error. Both route through surrogate_score() (R/likelihood.R), which divides max_batch down into cross_iid()'s per-block size; a NA/NaN max_batch (Inf was already fine) reached rep(..., times = NA) there and failed with the bare base-R message "invalid 'times' argument", naming neither max_batch nor log_lik()/log_ratio(). surrogate_score() now runs max_batch through check_positive(max_batch, "max_batch", allow_inf = TRUE) up front, matching the validation already applied to every other numeric argument at that boundary (#230) (#232).sample_prior() now validates n. It was the one count-taking public entry point that never routed its count through check_count(), unlike n_simulations, n_sbc, n_init, and simulate_for_sbi()'s own n. sample_prior(prior, 2.5) silently returned only 2 rows, via base R's recycling and an opaque warning; sample_prior(prior, -1) or sample_prior(prior, NA) failed with base R's generic "invalid arguments" error, naming neither the argument nor the function. sample_prior() now runs n through check_count(n, "n") at the top, before it reaches prior$sample() (#231).npe_sequential() now validates epsilon up front, before round 1 simulates. Every other round-controlling argument (n_rounds, n_simulations, n_truncation_samples, max_proposal_batches) is checked before any simulation runs, but epsilon fell through that net: it is only read starting in round 2, where it sets the truncation threshold via stats::quantile(lp_ref, probs = epsilon, ...). An out-of-range epsilon (e.g. 1.5, -0.1, NA) passed construction and round 1 silently, spent round 1's simulation budget, and only then failed in round 2 with stats::quantile()'s raw 'probs' outside [0,1] error. npe_sequential() now runs epsilon through check_prob() alongside the other up-front checks (#226) (#228).prior_uniform()/prior_normal() no longer let NA/Inf build a silently corrupted prior. Both constructors coerced their arguments with bare as.numeric() and never validated them, unlike every named family in R/prior_families.R. prior_uniform(low = c(mu = 0), high = c(mu = Inf)) built without error, and sample_prior() on it returned Inf on every draw, since runif(n) * (high - low) + low is always infinite once the range is; prior_uniform(low = NA, high = 5) crashed with a raw, unnamed base-R error ("missing value where TRUE/FALSE needed") instead of a validated one, because the high <= low comparison came before any check that low/high were finite; and prior_normal(mean = NA, sd = 1) built without error and silently returned an all-NA sample matrix. prior_uniform() now runs low/high through check_finite(allow_inf = TRUE), which rejects NA/NaN but still allows Inf -- needed so an improper prior can be built for prior_truncated()/prior_independent() to bound -- and prior_normal() runs mean/sd through check_finite() with no such allowance, since a normal prior has no legitimate use for an infinite mean or sd. Calling sample_prior() directly on an infinite-bound prior_uniform() now errors with a message pointing at prior_truncated(), rather than returning Inf (#227).c2st() no longer silently corrupts a whole column on a non-finite entry. x/y went through check_numeric() alone, and with the default z_score = TRUE a single NA/NaN/Inf anywhere in a column reaches fit_standardizer()/apply_standardizer() (R/standardize.R), whose colMeans()/sd() carry no na.rm -- so that column standardizes to NA in every row, not just the offending one. The corrupted matrix then either crashes several frames away with no mention of x or y (classifier = "logistic"'s glm() call errors with "Argument mu must be a nonempty numeric vector" once na.action drops every row; classifier = "mlp"'s training loop hits "missing value where TRUE/FALSE needed" once the propagated NaN reaches the stopping check), or, worse, survives to report an accuracy/AUC computed on a corrupted column, the metric this package's diagnostics are built around. c2st() now runs check_finite() on x and y right after check_numeric(), matching surrogate_score() (#202) and mcmc_log_prob() (#208/#209) (#222).log_prob() on an NPE posterior no longer returns a silent NaN for a non-finite theta. log_prob.nsbi_posterior() (R/posterior.R) validated theta with check_numeric() alone, never check_finite(). With a bounded prior, within_support() on a row containing NA/NaN returns logical NA, and R leaves an NA-indexed position untouched on assignment, so lp[!within_support(prior, theta)] <- -Inf was a no-op for that row and the NaN de_log_prob() produced reached the caller unnamed. log_prob.nsbi_posterior() now runs theta through check_finite(theta, "theta", allow_inf = TRUE) right after check_numeric(), matching mcmc_log_prob()'s precedent for the same "posterior log_prob" contract (#202) (#208) (#209): Inf stays allowed, since it resolves to zero density or -Inf through the prior or the estimator rather than signaling a bug, and only NA/NaN are rejected (#221) (#223).bulk_ess() no longer floors Geyer's tau estimate at 1. The reference formula (Stan's implementation, and Vehtari et al. 2021, which split_rhat() already matches) floors tau_hat at 1 / log10(n*k), not at 1; min(n*k / tau, n*k * log10(n*k)) and n*k / max(tau, 1/log10(n*k)) are the same clamp written two ways, and the code used the first form for the upper bound but then separately floored tau itself at 1, discarding the lower bound the reference formula allows. A well-mixing or anti-correlated chain routinely has an uncapped tau_hat below 1, and the wrong floor silently clamped its ESS at n*k, the raw draw count, instead of the larger value the reference computation gives -- confirmed against posterior::ess_bulk() on AR(1) chains with negative lag-1 autocorrelation, where the old code reported exactly n*k while the correct value ran 2x-3x higher. A new test adds anti-correlated AR(1) fixtures (#219) (#220).bulk_ess() no longer over-corrects autocorrelation at every lag. autocov() returns the biased autocovariance (divided by n), and the n / (n - 1) factor exists to convert only the lag-0 term into the unbiased quantity that matches W, which stats::var already computes with the n - 1 divisor. The code instead applied n / (n - 1) to the whole autocovariance vector before computing rho, which happened to leave lag 0 correct (acov[1] * n / (n - 1) == W algebraically) but inflated every higher-lag autocorrelation entering Geyer's paired sum, matching neither Stan's own implementation nor rstan's monitor.R::ess_rfun, which rescale acov[0] alone and leave every lag >= 1 unscaled. The effect shrinks as n -> Inf (the n / (n - 1) factor approaches 1), which is why the package's regression test against posterior::ess_bulk() didn't catch it: it ran at n = 400 with a 1% tolerance, small enough to hide a discrepancy that reaches 1-2% at the n = 20-100 range typical of nle()/nre()'s default slice-sampler settings. Only the diagnostic ESS number reported by mcmc_diagnostics() was affected, not the MCMC draws themselves or split_rhat(). A new test adds an n = 50 fixture with a tolerance tight enough to have caught this (#217) (#218).npe_sequential() now honors seed for weight initialization, not just for the proposal/simulation randomness. npe_sequential() calls set.seed(seed) up front, which seeds R's base RNG for the whole run, but its inner npe() call never forwarded seed on, so npe()'s own seed argument stayed at its default NULL. train_restarts() (R/train.R) only calls torch::torch_manual_seed(seed) when seed is non-NULL, and that call is what drives reproducible network weight initialization every round -- so two npe_sequential(..., seed = 42) calls could still train different networks, depending on whatever torch RNG state happened to be ambient at call time. The fix passes seed through to the inner npe() call (#216).npe()/nle()/nre() now honor seed when called with pre-computed theta/x. The seed argument only ever reached R's base RNG through simulate_for_sbi()'s own set.seed(seed), which runs on the branch of prepare_simulations() that calls the simulator; the documented, first-class theta/x code path skipped it entirely. train_restarts() (R/train.R) draws its train/validation split and every epoch's minibatch order from sample.int(), which reads R's base RNG, not torch's -- so two calls with an identical seed and identical pre-computed theta/x could still train on different splits and minibatch orders, and so land on different weights, depending on whatever base RNG state happened to be ambient at call time. All three entry points now call set.seed(seed) themselves before prepare_simulations() runs, matching the pattern already used in npe_sequential(), c2st(), and the NLE MCMC posterior (#214).de_sample.nsbi_de_mdn() no longer errors on a 1-d target with more than one mixture component. It read the batch row of the means and Cholesky-factor tensors with params$means[1, , ] / L[1, , , ], and R's torch package indexing follows base-R drop = TRUE semantics: every size-1 dimension in the result is dropped, not just the batch dim being integer-indexed. With dim_theta == 1 (or, under nle()'s swapped roles, dim_x == 1) and n_components > 1 -- the estimator's default is 10 -- the p dimension collapsed too, so means/Larr came back as bare vectors instead of K x p / K x p x p arrays and de_sample() errored with "incorrect number of dimensions". The existing n_components == 1L special case worked around the same drop for K, but nothing covered p == 1. The fix reshapes explicitly from the full, un-indexed tensor via array(torch::as_array(...), dim = c(K, p)) (and c(K, p, p) for Larr), which subsumes the old K == 1L branch, so it is removed (#211) (#212).c2st() now runs the same test as the sbibm benchmark, so its numbers are comparable with published ones. It fit a cross-validated logistic regression, which is linear: it scores two sample sets that share a mean and differ in spread at chance, and a difference in spread is exactly what a mis-trained posterior usually shows. R/c2st.R replaces it with the procedure in sbibm/metrics/c2st.py -- both sample sets z-scored by the mean and standard deviation of x, then a two-hidden-layer ReLU network of 10 * d units per layer trained by Adam (lr = 1e-3, L2 1e-4, minibatches of 200, stopping once the epoch loss has gone 10 epochs without improving by 1e-4), scored by accuracy over 5 shuffled folds. On two standard normals differing only in scale (sd 1 against sd 2) the old classifier returned 0.50 and the new one returns 0.74. The network trains on torch, like the rest of the package's neural code, so classifier = "mlp" needs torch installed and says so when it is missing; classifier = "logistic" is the torch-free way to get a number, and it is what the analytic-parity tests assert on in the CI job that has no libtorch.c2st() gains classifier, z_score, noise_scale, hidden, max_epochs and device. classifier = "logistic" is the old linear test, kept as a cheap screen. z_score and noise_scale are sbibm's arguments of the same name; the noise is for draws that are discrete or lie on a lower-dimensional set, where a classifier separates the two sides on an artefact of representation. device sends the network to a GPU the same way the fit_* functions do. It also returns the ROC AUC alongside the accuracy, as sbibm does.c2st()'s two sample sets are no longer symmetric, since x alone sets the z-scoring. Pass the reference draws as x, following sbibm. The call sites in tests/ and inst/benchmarks/ were flipped to match.inst/benchmarks/13_c2st_parity.R and 14_c2st_parity.py check c2st() against sbibm/metrics/c2st.py itself on five sample-set pairs whose answer is known by construction, so the parity claim rests on a run rather than on a reading of the Python source. The two agree to Monte-Carlo noise; inst/benchmarks/README.md records the table.log_prob() on an NLE/NRE posterior no longer returns a silent NA/NaN for a non-finite theta. mcmc_log_prob() (R/nle_posterior.R), the shared body behind log_prob.nsbi_mcmc_posterior(), validated theta with check_matrix() alone, which checks type and shape but never finiteness. An NA/NaN entry reached surrogate_potential()'s closure, came back non-finite from prior$log_prob() the same way an out-of-support draw does, and was left as a silent NA/NaN log-prob instead of a named error. mcmc_log_prob() now runs theta through check_numeric() and check_finite() before check_matrix(), matching surrogate_score() (#202) and stan_data() (#203). check_finite() gains an allow_inf argument for this call site: an Inf theta is not part of the bug, since it correctly resolves to -Inf through the prior's own density, so it stays allowed and only NA/NaN are rejected (#208) (#209).stan_data() now validates x_obs before handing it to Stan. posterior.nsbi_npe() and mcmc_posterior() both check x_obs with check_numeric() and check_finite() before as_theta_matrix() (#51), but stan_data() called as_theta_matrix() directly. An NA in an otherwise numeric x_obs, or a character column that storage.mode<- silently turns into NA, reached cmdstan_model()$sample()/rstan::sampling() unnoticed and surfaced as an opaque Stan error instead of a named neuralsbi one. stan_data() now runs the same two checks posterior() and mcmc_posterior() already do (#203) (#206).log_lik() and log_ratio() no longer return a silent NA for a non-finite theta or x. surrogate_score() (R/likelihood.R), the body both share, validated its arguments with check_matrix() alone, which checks type and shape but never finiteness. An NA/NaN/Inf entry standardized into another NA and came back as a silent NA log-density or log-ratio instead of an error, the same failure mode fixed for theta passed to a surrogate posterior's log_prob() in #163. surrogate_score() now runs theta and x through check_numeric() and check_finite() before check_matrix(), matching posterior.nsbi_npe() and mcmc_posterior() (#202) (#205).
stan_data() now validates x_obs before handing it to Stan. posterior.nsbi_npe() and mcmc_posterior() both check x_obs with check_numeric() and check_finite() before as_theta_matrix() (#51), but stan_data() called as_theta_matrix() directly. An NA in an otherwise numeric x_obs, or a character column that storage.mode<- silently turns into NA, reached cmdstan_model()$sample()/rstan::sampling() unnoticed and surfaced as an opaque Stan error instead of a named neuralsbi one. stan_data() now runs the same two checks posterior() and mcmc_posterior() already do (#203) (#206).
tests/testthat/test-stan.R checks that the Stan code stan_code() emits for linear_gaussian, mdn, and maf agrees numerically with log_lik(), but skip_if_no_cmdstan() (tests/testthat/helper-stan.R) meant they ran nowhere: neither R-CMD-check.yaml nor test-coverage.yaml installs a CmdStan toolchain. A new .github/workflows/stan-tests.yaml installs libtorch and CmdStan (both cached across runs) and runs the Stan tests on workflow_dispatch and on a schedule every two weeks, on main, rather than on every push -- installing a compiler toolchain on the fast-feedback workflows would slow down every push and PR for a job that only exercises one file (#196) (#204).prior_lognormal(), prior_exponential(), prior_gamma(), prior_beta(), prior_student_t(), prior_cauchy(), prior_half_normal() and prior_half_cauchy() (R/prior_families.R) build a prior from a named distribution family under Stan's argument names, vectorized over parameters the way prior_normal() already is. A named family sets lower/upper from its own support, so the out-of-support rejection and the log_prob renormalization in R/posterior.R get the right region with nothing further to declare (#199) (#201).prior_independent() multiplies per-parameter priors into a joint one. This is what most prior_custom() calls were written to do by hand, and doing it by hand is where the quiet mistakes live: a log_prob_fn that returns one number instead of one per row, or a lower of the wrong length that sweep() recycles into a support test rejecting the wrong draws. Components may be any nsbi_prior, including multi-parameter ones and a prior_custom(). Argument names name one-parameter components; a wider component keeps its own param_names (#199) (#201).prior_truncated() is Stan's T[lower, upper], with the density renormalized by the mass it keeps. Leaving the normalizing constant off would be harmless for a posterior sampled by MCMC on its own, but nle() and nre() sum the prior density with a learned likelihood in surrogate_potential() (R/mcmc.R), so a prior short a constant is a prior of the wrong shape relative to that likelihood. Bounds intersect with the family's own support and with any earlier truncation. Truncating a prior_custom() errors rather than returning an unnormalized density, since there is no CDF behind it to renormalize against (#199) (#201).stan_code() writes the new families out as sampling statements. Uniform and normal priors still travel through the data block, unchanged. Every other family becomes a literal sampling statement carrying T[,] on whichever side the support was cut, and the parameter block declares matching constraints: a shared bound gives vector<lower=...>[Q] theta, differing bounds give per-parameter reals assembled in transformed parameters. A prior_custom() still errors, now naming ?prior_families as the alternative (#199) (#201).prior_uniform() and prior_normal() now carry the same internal per-parameter form as the new families, so both can be truncated and composed. Their closures, printing and Stan output are unchanged. print.nsbi_prior() gains one line per marginal for the new types (p ~ beta(2, 15), s ~ normal(0, 1.5) T[0, ]) (#199) (#201).task_sir()'s prior (R/tasks.R) is now a prior_lognormal() rather than the same two log-normals written out through prior_custom(), so an nle() fit on that task exports with a model block instead of erroring (#199) (#201).vignette("methods"), a guide to the methods this package implements: NPE, sequential NPE, NLE and NRE, the MDN/MAF/NSF density estimators, and the SBC, coverage, TARP and C2ST diagnostics. Each entry names the paper the method comes from and then the function that runs it, so a reader who knows a method by its citation can find the call. Citations resolve through inst/REFERENCES.bib, so every entry renders with a link out to the DOI (#200) (#201).mcmc_init(strategy = "resample") no longer double-counts the prior in its SIR weight. mcmc_init_resample() (R/mcmc.R) draws its candidate pool from the prior -- so the SIR proposal is the prior itself -- and resampled without replacement using log_prob_fn(cand) as the weight, which for nle()/nre() posteriors is surrogate_potential()'s full unnormalized posterior, log p(theta) + log L(x|theta). A correct importance weight is target/proposal, and since target is p(theta) * L(x|theta) and the proposal is p(theta), the prior cancels and the weight should be L(x|theta) alone; weighting by the full posterior instead resampled from p(theta)^2 * L(x|theta), pulling chain starts back toward the prior's mode. prior_uniform() hid this, since a constant drops out of the Gumbel-top-k ranking either way -- the bug only shows for a non-flat prior, such as prior_normal() or a prior_custom() with an asymmetric density. mcmc_init_resample() now subtracts prior$log_prob(found) from the accumulated found_lp before building the Gumbel keys, leaving the finiteness check that selects found on the full posterior density unchanged (#195) (#197).npe_sequential() now errors when n_simulations doesn't match n_rounds, instead of silently recycling it. n_simulations is documented as "either a scalar or a vector of length n_rounds" (R/sequential.R), but check_counts() only validated each element, never length(n_simulations) against n_rounds; budgets <- rep_len(n_simulations, n_rounds) then recycled any other length without a diagnostic. npe_sequential(prior, sim, n_rounds = 3, n_simulations = c(500, 1000)) -- a plausible mistake, such as forgetting the third round's budget -- silently ran with budgets c(500, 1000, 500) instead of raising an error. Since every round retrains from scratch and the point of the multi-round scheme is spending a fixed budget deliberately across rounds, a silently wrong per-round budget wastes or misallocates simulator calls with nothing to catch it. npe_sequential() now checks length(n_simulations) %in% c(1L, n_rounds) right after check_counts() validates its elements, and errors naming both arguments before rep_len() runs (#191) (#194).max_steps and n_pool, closing the gap the width fix in 0.5.26 explicitly left open. slice_sample_run() (R/mcmc.R) fed a user-supplied max_steps straight into steps_left <- max_steps; while (steps_left > 0L), and mcmc_init()'s n_pool fed straight into batch <- max(n_pool, n_chains) inside mcmc_init_resample()/mcmc_init_proposal(), neither checked before use. Both are reachable from posterior(nle_fit_or_nre_fit, x_obs, max_steps = ..., n_pool = ...): posterior.nsbi_nle()/posterior.nsbi_nre() forward them through ... to slice_sample_surrogate() (R/nle_posterior.R), which passes dots$max_steps %||% 100L to slice_sample() and dots$n_pool %||% 1000L to mcmc_init(). max_steps = NA looped forever instead of erroring, since NA > 0L is NA and a while() condition only errors once evaluated; a zero or negative n_pool reached max() and sample_prior() with nothing said about which argument was wrong. slice_sample_run() and mcmc_init() now check both with check_count() (R/check.R), right where check_slice_width() already checks width, before either value can drive a loop or size a batch (#190) (#192).nre() no longer silently keeps an untrained network when a small n_simulations leaves a 1-row validation split. check_train_controls() (R/train.R) only required the validation split to be non-empty, which is enough for npe()/nle() -- their estimators score a real, if noisy, log-density on a single validation row. nre()'s atomic contrastive objective (nre_atomic_log_prob(), R/nre.R) cannot: with fewer than 2 rows there is no contrast to score, and it returned a literal constant-zero loss every epoch. Inside train_restarts() that zero beat the initial Inf at epoch 1, never changed afterward, so the "improved" branch never fired again, and training silently stopped after patience more epochs holding the epoch-1 (essentially untrained) network -- while reporting best_val_loss = 0, which reads as a perfect fit rather than the failure it is. nre(prior, simulator, n_simulations = 15) hit this directly: the default validation_fraction = 0.1 gives n_val = max(1, floor(0.1 * 15)) = 1. check_train_controls() now takes a min_val_rows argument (default 1L, unchanged for every other caller) threaded through train_conditional_de() and fit_torch_de(); fit_nre_net() (R/nre.R) passes min_val_rows = 2L, and nre() checks it early, before the simulator runs, whenever the true row count is already known (theta/x passed directly, or a valid n_simulations). The closed-form logistic classifier and a caller-supplied classifier function are unaffected -- neither uses this validation split at all (#188) (#189).map_estimate() no longer optimizes a 1-parameter posterior with Nelder-Mead. stats::optim()'s own documentation and runtime warn that one-dimensional Nelder-Mead is unreliable ("use "Brent" or optimize() directly"), and single-parameter models are not a corner case for this package -- one-parameter priors appear in the nle()/nre() doc examples themselves. Every map_estimate() call on such a posterior hit that warning, and a test already worked around it with suppressWarnings() instead of fixing the underlying method choice (tests/testthat/test-posterior-normalization.R). map_estimate() (R/posterior.R) now branches on fit$dim_theta == 1L: a bounded prior with both a lower and an upper limit gets optim(method = "Brent") over that exact interval; a one-sided bound gets optim(method = "L-BFGS-B") with Inf on the missing side, since Brent needs both ends finite and plain BFGS's finite-difference gradient can probe past the missing side into masked territory; a fully unbounded prior gets optim(method = "BFGS"), which needs no interval and doesn't trigger the warning. Multi-dimensional posteriors are unaffected and stay on Nelder-Mead (#186) (#187).mcmc_init(strategy = "resample") no longer starts several chains from identical points. It drew one pool of max(n_pool, n_chains) prior draws and kept whichever landed inside the posterior's support; when that left fewer than n_chains of them, it padded the shortfall with rep_len(ok, n_chains), recycling the same indices so several slice-sampler chains began from the same point (R/mcmc.R). Split-Rhat and bulk ESS both assume the chains they compare started from distinct locations, so a recycled start weakens mode coverage and understates disagreement between chains. A new mcmc_init_resample() draws more pools and accumulates their finite draws, the way mcmc_init_proposal() already does for "proposal", until n_chains distinct finite draws are collected (up to 20 attempts) and only then runs the weighted resample without replacement. If the budget runs out first, the error distinguishes a genuinely degenerate surrogate (no draw at all landed in support) from a merely low acceptance rate (some did, just not enough) (#182) (#185).standardizer_log_jac()'s roxygen docstring (R/standardize.R) named the wrong transform direction. It said "inverse standardization (standardized -> original)", but the value it returns, -sum(log(std$scale)), is the Jacobian of the forward transform z = (x - center) / scale (original -> standardized); the code was always correct, every call site already adds it to a standardized-space log-density to recover original units. The docstring and man/standardizer_log_jac.Rd now describe the forward direction (#183) (#184).mcmc_init(strategy = "proposal") no longer fails on an ordinary, non-degenerate posterior with a modest acceptance rate. It drew one batch of n_chains prior draws and required all of them to be jointly finite, retried whole-batch up to 20 times (R/mcmc.R). If the posterior excludes even a third of the prior's support -- an entirely ordinary case, not a degenerate one -- a batch is finite with probability 0.65^20, about 3e-4, so the retries almost never succeeded, and the error claimed the surrogate likelihood "may be degenerate" regardless of why the batch failed. A new mcmc_init_proposal() pools finite draws across attempts instead, the way "resample" already pools its candidates, so the chance of collecting n_chains finite starting points grows with the total number of prior draws tried rather than with the luck of one batch. If the budget still comes up short, the error now distinguishes a truly degenerate surrogate (no draw at all landed in support) from a merely low acceptance rate (some did, just not enough), and only the former says "degenerate" (#179) (#181).npe_sequential() no longer crashes with an opaque rbind() error when a later round accepts zero proposals and dim_x > 1. The truncated-proposal loop (R/sequential.R) calls run_simulator() for round r's accepted theta_new, but didn't pass its optional d argument (expected output width). When a round exhausted max_proposal_batches with zero accepted candidates, run_simulator()'s zero-row fast path (R/parallel.R) fell back to d %||% 1L and returned a 0 x 1 matrix regardless of the simulator's real output width, so x_all <- rbind(x_all, x_new) failed with "number of columns of matrices must match" whenever dim_x != 1 -- a message naming neither npe_sequential() nor the real cause. With dim_x == 1 it didn't crash, but silently contributed zero simulations to that round with only a warning to go on. npe_sequential() now passes d = ncol(x_all) once round 1 has run, and treats zero accepted proposals as a hard stop rather than "continuing with fewer simulations": a round that adds nothing would otherwise refit on unchanged data and report that as round r's result, with no way for a caller to see that nothing happened. The new error names the round, the acceptance, and suggests a larger epsilon, more max_proposal_batches, or fewer rounds (#178) (#180).expected_coverage() now returns the right shape for a single-parameter fit. For a fit with one parameter, sbc_result$ranks has one column, so colMeans() inside expected_coverage()'s (R/diagnostics.R) per-level sapply() returned a length-1 result at every nominal level; sapply() then simplified those to a plain length(levels) vector instead of a matrix, and the following t() turned that into a single-row matrix -- one column per nominal level instead of one column for the parameter. colnames() then labeled those columns param1/param2/param3 as if there were three parameters, each holding the same three numbers repeated across every nominal level, so the table read as flat (and wrongly shaped) regardless of how well calibrated the fit actually was. expected_coverage() now uses vapply() and restores the dimension it drops for a one-column ranks matrix before transposing, matching the shape it already produced correctly for two or more parameters. Found while investigating #169: a one-parameter fit is not what that issue's SIR model uses (three parameters), so this bug is not the cause of the undercoverage reported there, but it is a real, previously uncaught defect in any single-parameter calibration check (#169) (#177).sbc()'s calibration under a uniform prior that actively truncates the fitted density -- a linear_gaussian fit with its B/Sigma set by hand to the exact conditional (theta | x ~ N(x, sigma^2 I) restricted to the prior box, the closed form for a flat prior and additive Gaussian noise), so the only code path left to trust is the rejection-sampling and acceptance-renormalization in R/posterior.R and the rank binning in sbc() itself. No existing test checked SBC calibration with an exact fit under an actively truncating bounded prior: the existing calibration tests all use an unbounded prior_normal(), and the existing bounded-prior tests only use prior_uniform() to check the "short draw errors" behavior with a fit deliberately broken to leak. This test comes back well calibrated, which rules out both of #169's bounded-prior-handling and rank-scoring candidates as the source of that issue's undercoverage (#169) (#177).posterior(fit, sampler = "stan") now falls back to rstan (or says clearly what to run) instead of crashing when cmdstanr has no CmdStan behind it. stan_sample_nle() (R/stan.R) picked its backend with a bare requireNamespace("cmdstanr", quietly = TRUE), which only checks that the R package is installed. install.packages("cmdstanr") does not install CmdStan itself -- that is the separate, network-dependent cmdstanr::install_cmdstan() step, which fails outright on machines that block the GitHub release download it needs (many CI runners and managed environments). In that state the old code always took the cmdstanr branch and failed deep inside cmdstanr::cmdstan_model() with "CmdStan path has not been set yet.", even when a working rstan install was sitting right there. A new cmdstan_ready() helper checks cmdstanr::cmdstan_version(error_on_NA = FALSE) as well as the package, mirroring the package-vs-runtime distinction require_torch() already draws for the torch backend. stan_sample_nle() now uses it to pick the backend: it falls back to rstan (with a message explaining why) when cmdstanr is installed but CmdStan is not, and only errors when there is truly nothing to fall back to -- with a message that tells cmdstanr-without-CmdStan apart from neither package being installed, since the two have different fixes (#172) (#176).posterior_predictive() now errors clearly when the posterior has no draws inside the prior support, instead of crashing on a dimnames assignment. For a bounded prior, sample.nsbi_posterior() can legitimately return zero rows when rejection sampling never lands inside the support for an observation far outside anything the fit was trained on -- it warns, but does not stop. posterior_predictive() (R/diagnostics.R) passed that zero-row theta straight to run_simulator(), which produced a 0x0 matrix, and the function then ran colnames(pred) <- post$fit$x_names against it: with dim_x > 1 that failed with the generic length of 'dimnames' [2] not equal to array extent, naming neither the function nor the actual problem, and with dim_x == 1 it silently returned a useless empty matrix instead of erroring at all. posterior_predictive() now checks nrow(theta) right after sampling and stops with a message naming the function, the acceptance rate, and the likely cause (the observation being outside the range the fit was trained on) before the simulator or the dimnames assignment ever run. summary.nsbi_samples() (R/summaries.R) gets the same treatment for the related silent failure the issue flagged: summarizing a zero-row draw returned a row of NaN/NA per parameter with no warning, and now says so (#171) (#174).width instead of letting a bad value crash the stepping-out loop several frames later. nle()/nre() posteriors forward a width argument through ... from posterior() down to slice_sample()/slice_sample_run() (R/mcmc.R), which coerces it with width <- rep_len(as.numeric(width), dim) and never checked the result before using it to size the initial slice interval and drive the stepping-out loop's array indexing. posterior(nle_fit, x_obs, width = NA) followed by sample() failed with "NAs are not allowed in subscripted assignments" -- a generic R subscripting complaint naming neither width nor the sampler. slice_sample_run() now checks the recycled width with a new check_slice_width() (R/mcmc.R) right after rep_len(), and stops with a message naming width, the values that were actually wrong, and slice_sample(), before any of them reach stats::runif() or the stepping-out loop. Negative/zero width and negative max_steps/zero n_pool degrade silently rather than crash and are not covered by this fix (#164) (#167).log_prob() on an NLE or NRE posterior now rejects a non-numeric theta by name instead of silently returning NA. log_prob.nsbi_mcmc_posterior() (R/nle_posterior.R) routes theta through mcmc_log_prob() into surrogate_potential()'s closure (R/likelihood.R), which only reshapes it with as_theta_matrix() -- no check_matrix()/check_numeric() anywhere on that path. as_theta_matrix() coerces a non-numeric column to NA with a bare "NAs introduced by coercion" warning rather than erroring, so log_prob(post, theta = data.frame(mu = c("x", "y"))) came back as NA with nothing said about theta -- worse than crashing, since there is no signal to a caller doing which.max(log_prob(...)). Every sibling entry point already validates: log_prob.nsbi_posterior() for NPE posteriors, and log_lik()/log_ratio() for the same NLE/NRE fits via surrogate_score(). mcmc_log_prob() now checks theta with the same check_matrix() call surrogate_score() already uses, once, before building the potential -- not inside surrogate_potential()'s returned closure, which every MCMC step also calls with a theta it built internally and has no need to re-validate (#163) (#166).sample() now validates n (aliased size), matching every other draw-count argument in the package. sample.nsbi_posterior() (R/posterior.R) already checked max_sampling_batches with check_count(), but n -- the argument that actually drives the rejection-sampling loop -- reached while (nrow(collected) < n ...) unchecked, so sample(post, n = NA) failed inside the loop condition, sample(post, n = "10") failed on n - nrow(collected), and sample(post, n = -5) failed inside matrix(), none of the messages naming n or sample(). sample.nsbi_mcmc_posterior() (R/nle_posterior.R, backing nle()/nre() posteriors) had the same gap and it was worse there: n reached mcmc_draws()'s ceiling(n_draws / n_chains) unchecked, so sample(post, n = 2.5) did not error at all, it silently returned 2 draws. Both methods now check n with check_count() up front, before any sampling work happens, giving the same named, actionable error n_init, n_sbc and the rest of the package's draw-count arguments already give. posterior_predictive() (R/diagnostics.R), which forwards its own n straight into sample(), inherits the fix without any change of its own (#162) (#165).map_estimate() now errors clearly instead of crashing when its seeding draw comes up short. It seeds stats::optim() with draws <- sample(post, n = n_init, obs = x) and then picks the best of them by log_prob(), but sample.nsbi_posterior() can legitimately return fewer than n_init rows -- including zero -- for a bounded prior when rejection sampling never lands inside the support after max_sampling_batches rounds, and it only warns, it does not stop. A zero-row draw reached which.max() and then stats::optim() unchecked, and the failure surfaced deep inside dmvnorm_chol() (R/density_estimator.R) as Error in x - mean : non-conformable arrays, naming neither map_estimate() nor the leaking prior actually responsible. map_estimate() (R/posterior.R) now checks nrow(draws) right after the sample() call, matching the check diagnostic_draws() (R/diagnostics.R) already makes for the same failure mode, and stops with a message naming the function and the shortfall and suggesting more simulations or a look at the prior. A partial short draw errors too, not just an empty one -- there is no way to know n_init points were searched from when fewer were (#159) (#161).split_rhat() now includes the tail-Rhat component, so it matches its Vehtari et al. Rhat docstring in full. #157 fixed the missing rank-normalization but stopped at bulk-Rhat: it rank-normalized the raw split-chain draws and ran the classical Gelman-Rubin statistic on those, without also folding them around the median for a tail component. Vehtari et al. (2021)'s Rhat -- and posterior::rhat()'s default -- is max(bulk-Rhat, tail-Rhat); bulk-Rhat alone tracks location but not scale, so it reads chains that agree in location but disagree in spread as converged. Four chains centered at 0 with two at sd = 1 and two at sd = 4 reproduce this: bulk-Rhat alone reads Rhat < 1.01 while the complete statistic reads Rhat > 1.2. split_rhat() (R/mcmc.R) now computes both components through a shared gelman_rubin_rhat() helper and takes their max, and tests/testthat/test-mcmc.R's posterior::rhat() cross-check tightens from tolerance = 1e-4 back to 1e-6 now that the two compute the same statistic rather than a close approximation of it (#154) (#160).split_rhat() now actually rank-normalizes, matching what its docstring already claimed. mcmc_diagnostics() documented "the standard rank-free versions from Vehtari et al. (2021)" for both statistics it returns, and bulk_ess() did rank-normalize, but split_rhat() computed the classical 1992 Gelman-Rubin statistic straight off the raw split-chain values. Rank-normalization is what makes Rhat robust to heavy tails and multimodality -- exactly the failure mode the vectorized slice sampler (R/mcmc.R) behind nle()/nre() posteriors is prone to -- so a raw-value Rhat can read as converged when the rank-normalized one would flag it: a four-chain, Cauchy-noised fixture with real between-chain disagreement reads Rhat < 1.01 under the old code and Rhat > 1.2 under the fix. split_rhat() and bulk_ess() now share one rank_normalize() helper, so there is one implementation of Vehtari et al.'s rankit transform instead of the claim living in a docstring the code did not follow. tests/testthat/test-mcmc.R gains the Cauchy regression fixture and tightens its posterior::rhat()/posterior::ess_bulk() comparison tolerance now that the two are computing the same thing rather than merely landing close on near-Gaussian data (#154) (#157).log_prob(..., normalize = TRUE) now warns when its acceptance estimate hits the zero floor, instead of staying silent. For a bounded prior, log_prob.nsbi_posterior() (R/posterior.R) estimates the normalizing acceptance probability by sampling n_normalization draws and checking within_support(); when none of them land inside the prior, it floors the estimate at 1 / n_normalization to avoid log(0) and returns a normal-looking finite value with no sign that the normalizing constant is a floor rather than a measurement. sample.nsbi_posterior() already warns in the equivalent case, when rejection sampling comes up short: "The estimator is leaking mass outside the prior; consider more simulations." log_prob() now raises the same warning whenever the floor is doing the work, so the two functions agree about surfacing this failure mode (#153) (#156).map_estimate() no longer returns a point outside a bounded prior's support. It searches with unconstrained stats::optim(..., method = "Nelder-Mead") and always calls log_prob(post, ..., normalize = FALSE), so the -Inf mask that log_prob(normalize = TRUE) applies for a bounded prior (from prior_uniform() or a prior_custom() with lower/upper) never reached the optimizer's objective, and the search could converge outside the box. The objective in map_estimate() (R/posterior.R) now checks within_support() directly and returns Inf for an out-of-support candidate, which is the same guarantee without paying to re-estimate the acceptance constant on every evaluation. Unbounded priors take the same code path as before and are unaffected (#152) (#155).#136 fixed the parts of #107 that don't need libtorch to verify (pkgdown's equation rendering and pairplot()'s truth-marker/axis-range styling); this closes the remaining part, converting every base-plot() call across vignettes/*.Rmd.orig (intro-to-sbi, neural-likelihood, neuralsbi, sir-epidemic, sir-time-varying-beta) to ggplot2, matching the theme_minimal() and steelblue/firebrick/grey palette the package's own plot_*() functions already use. neural-likelihood.Rmd.orig's SBC panel loop also drops a par(mfrow = c(2, 2)) wrapper around four calls to plot_sbc(), which is grid graphics and was never affected by par(); the baked figures already rendered as four separate images, so removing the dead wrapper changes nothing about the output. These vignettes are precomputed under libtorch (see vignettes/precompute.R) and this environment has no libtorch, so the committed figures under vignettes/figures/ are not yet regenerated from the new code -- that rebake is tracked in #29 (#107) (#148).sample(), log_prob() and print() each had a byte-identical method for nsbi_nle_posterior and one for nsbi_nre_posterior, differing only in the class name they printed and, for log_prob(), the "NLE"/"NRE" label in its unnormalized-posterior warning. mcmc_posterior() (R/nle_posterior.R) now stamps a shared nsbi_mcmc_posterior class into the middle of the class vector -- ahead of nsbi_posterior, behind the specific nsbi_nle_posterior/nsbi_nre_posterior that print() output and expect_s3_class() checks still key off -- so the three pairs collapse into sample.nsbi_mcmc_posterior(), log_prob.nsbi_mcmc_posterior() and print.nsbi_mcmc_posterior(). log_prob.nsbi_mcmc_posterior() reads the "NLE"/"NRE" label off inherits(post$fit, "nsbi_nle") instead of taking it from the caller, and print.nsbi_mcmc_posterior() gates its closing stan_code() line the same way. posterior.nsbi_nle() and posterior.nsbi_nre() are unchanged (#144) (#147).nre() adds Neural Ratio Estimation, the third factorization of the joint. Where npe() learns the posterior and nle() learns the likelihood, nre() (R/nre.R) learns neither density but their ratio r(theta, x) = p(x | theta) / p(x), by training a classifier to tell (theta, x) pairs drawn from the joint apart from pairs whose parameter came from a different simulation. Defaults mirror Python sbi's NRE (an alias for its NRE_B): the "resnet" classifier (a residual MLP, hidden = 50, n_blocks = 2), the atomic loss of Durkan et al. (2020) at num_atoms = 10, and the training controls npe()/nle() already share with sbi (batch_size = 200, lr = 5e-4, validation_fraction = 0.1, patience = 20, clip_grad_norm = 5). "mlp" and "linear" are sbi's other two classifiers; "logistic" is a torch-free baseline. log_ratio(fit, theta, x) evaluates the learned ratio, summing over rows of x as independent observations of the same parameter the way log_lik() does for an NLE fit, and posterior(fit, x_obs) returns an nsbi_nre_posterior that samples r(theta, x) p(theta) with the same vectorized slice sampler an NLE posterior uses. sbc(), tarp(), expected_coverage(), posterior_predictive(), map_estimate(), summary() and save_nre()/load_nre() all take the new fit. There is no stan_code() for a ratio estimator: the Stan exporter transpiles a fitted density, and a residual classifier has no such export."logistic" classifier is to nre() what "linear_gaussian" is to npe()/nle(): a closed-form, torch-free baseline that is exact for a linear-Gaussian simulator and so serves as the CI oracle for the whole NRE path. It fits a ridge-penalized logistic regression on the quadratic basis [1, z, vech(z z')] for z = (theta, x). Two atoms make the atomic loss collapse to a logistic regression on the difference of two feature rows, which is what IRLS can solve in closed form; working in differences also cancels every term depending on x alone, including the evidence log p(x), which is the one part of the log ratio a quadratic basis cannot represent. Contrasts are cyclic shifts rather than random draws, so the fit is deterministic given the simulations.p(x | theta) and p(x), so log_ratio() in the estimator's z-scored space is already the ratio in the units the simulator returned. surrogate_ops() (R/likelihood.R) is the one table where each surrogate fit type says which functions score it and whether the Jacobian travels with them, so surrogate_potential() (which replaces the internal nle_potential()), log_lik() and log_ratio() all read the answer off it instead of each restating it.linear_gaussian's ridge no longer depends on the scale of the data. fit_linear_gaussian() (R/density_estimator.R) added a flat 1e-06 to the diagonals of X'X and of the residual covariance. That is only a negligible regularizer when the data are O(1), which is what standardize = TRUE arranges and why the default path was never affected. Under standardize = FALSE, a target column whose variance is 2.5e-07 had 1e-06 added to it and came back five times too wide, throwing the reported log-density off by 20 nats on a three-column example spanning 1e-03 to 1e+03. The ridge is now measured per column against that column's own scale: diag(X'X) for the design, and the target's variance for the covariance, since a noiseless model leaves no residual to measure against and that is exactly when the ridge has to keep chol() alive. A column with no scale of its own borrows the largest one that has any. On standardized data both diagonals are already about 1, so the fits every other test pins are unchanged to six significant figures (#139).nre() fit. One simulation has no contrasting parameter to be scored against, so the atomic loss returned a constant that torch had not built and backward() failed with "element 0 of tensors does not require grad and does not have a grad_fn". Whether it happens at all depends on nothing but n_simulations modulo batch_size: 2001 simulations leave 1801 training rows, which at the default batch_size = 200 is nine full batches and one row. minibatches() (R/train.R) now folds a trailing single row into the batch before it, for every estimator rather than just the ratio ones -- a gradient step from one sample is a bad step for any of them -- and the atomic loss returns a tensor torch built, so a validation split of one row cannot reintroduce the failure (#139)."logistic" classifier's ridge no longer depends on the scale of the data either. irls_logistic() (R/nre.R) added a flat 1e-06 to the diagonal of X'WX, and its features are quadratic in (theta, x), so under standardize = FALSE they carry the fourth power of the data's units. On a simulator whose output has sd 5e-04 that ridge swamped the normal equations and the fitted ratio collapsed to noise: RMSE 5.88 against an analytic log ratio whose own spread is 6.29. The ridge is now measured per column against that column's own scale, through the same ridge_scale() that fit_linear_gaussian() uses, and the standardize = FALSE fit matches the standardize = TRUE one to six decimal places (#139).iid_matrix()/iid_evaluator() blocking loop is now tested without torch. linear_gaussian has its own de_log_lik_iid() and de_iid_evaluator() methods that score the whole observation set per parameter and ignore max_batch, so the existing "blocking does not change the answer" test compared two calls that did the same thing; cross_iid(), the path every flow and every ratio estimator takes, was reached only by the torch tests that CI skips. A scorer that encodes which (theta, x) pair it was handed now pins the block layout across eight batch sizes, so a transposition cannot hide inside a row sum (#139)."resnet" classifier is now pinned to nflows' ResidualNet numerically rather than by comment. sbi's "resnet" is that network, and nre_module() reimplements it; filling both with the same deterministic weights and evaluating them on the same rows agrees to 7e-08, which is float32. test-nre.R carries the three reference numbers, so where the activations sit inside a block, which of its two linear layers is zero-initialized, and whether the final layer sees an activation are all checked against the original (#139).nle() and nre() posteriors share are now written once. log_lik() and log_ratio() had the same body apart from the Jacobian and the pair of functions that score the estimator, so both now call surrogate_score() (R/likelihood.R), which reads them off surrogate_ops(). npe(), nle() and nre() built the same thirteen-field list in the same order, and now call new_nsbi_fit() (R/npe.R) with whichever fields are their own. cross_iid() (R/likelihood.R) takes the per-pair scorer as an argument, so de_log_lik_iid()'s and de_iid_evaluator()'s default methods and their ratio counterparts split into iid_matrix()/iid_evaluator() over one blocking loop; R/nle_posterior.R grows mcmc_posterior(), mcmc_draws(), mcmc_log_prob() and cat_mcmc_posterior() for the argument checking, construction, sampling, log-density and printing that both posterior classes do identically, and slice_sample_nle() becomes slice_sample_surrogate(); summary.nsbi_npe()'s body moves to fit_summary(), which the NRE method calls with classifier where the others pass density_estimator.sample.nsbi_posterior()'s support-rejection loop (R/posterior.R) and npe_sequential()'s truncated-proposal loop (R/sequential.R) both looped while (collected < n) { draw n more; keep the accepted ones }, asking de_sample()/sample_prior() for a full n (or full round budget) every round regardless of how many draws were already collected. For an estimator with real leakage at, say, 60% acceptance, requesting size = 1000 cost roughly 1000 + 1000 + ... per round instead of 1000 + 400 + ..., wasting sampling and support-checking work on later rounds. Both loops now request n - nrow(collected) (or budgets[r] - nrow(theta_new)) each round instead of the full count. Neither loop's stopping condition or returned shape changes -- both already stopped once n accepted draws were collected -- so this is a pure efficiency fix, not a behavior change (#138) (#143).npe()/nle() reject n_bins < 2 for the NSF density estimator instead of crashing on the first forward pass. nsf_made_module()'s forward splits a per-dimension parameter tensor into K bin widths, K bin heights and K - 1 interior derivatives by slicing params[, , (2 * K + 1):(3 * K - 1)]; for K = 1 that range is 3:2, which R's : reads as the descending index c(3, 2) rather than an empty selection, and the parameter tensor only has 2 columns for K = 1, so the read went out of bounds. check_architecture() already validated n_bins with check_count(), but its default min = 1L let n_bins = 1 through; it now passes min = 2L, the honest lower bound for a spline that needs an interior derivative to fit. Closes #137 (#140).\eqn/\deqn markup in man/*.Rd reaches reference/*.html through a different pkgdown code path than vignettes and README do: reference pages are built straight from Rd XML (as_html.tag_eqn/as_html.tag_deqn) and never go through pandoc, so the LaTeX landed on the page as literal \(...\)/$$...$$ text with nothing to render it, while vignettes (built from actual Markdown, through pandoc) rendered correctly under pkgdown's mathml default. _pkgdown.yml now sets template: math-rendering: katex, which works client-side by scanning the whole rendered page for math delimiters regardless of which pipeline produced the surrounding HTML -- fixing reference pages without changing how vignettes already render (#107) (#136).pairplot()'s truth markers are grey and dashed instead of a solid saturated red, and its panels now share a consistent axis range per parameter. The reference-value lines and cross marker used colour = "firebrick"; they're now "grey30", thinner, and dashed, so they read as an overlay rather than competing with the density fill. Separately, when limits isn't supplied, lower_fn/diag_fn each built an independent ggplot() object per panel, and ggdensity::geom_hdr()/ggplot2::geom_density() estimate their density grid per panel, so a parameter's drawn range could disagree between its diagonal panel and its off-diagonal panels, or between panels in the same row/column. pairplot() now defaults limits to each parameter's own data range (padded 5%, matching ggplot2's own default scale expansion) and applies it through the same coord_cartesian() path the explicit-limits argument already used, so every panel for a given parameter agrees (#107) (#136).sbc() no longer emits "Chi-squared approximation may be incorrect" during ordinary use. Its per-parameter uniformity test called stats::chisq.test(tab) on the binned ranks, and R warns whenever a cell's expected count is low, which happens routinely at the small n_sbc a test suite uses to keep runtime down. Switching to chisq.test(tab, simulate.p.value = TRUE) computes a Monte Carlo p-value instead of relying on the asymptotic approximation, so the warning never fires; the p-value keeps the same "large = calibrated" interpretation expected_coverage() and existing tests rely on (#109) (#135).R/utils.R, R/progress.R, R/simulator.R, R/parallel.R, R/mcmc.R, R/sequential.R, R/tasks.R, R/diagnostics.R and R/check.R -- mocked requireNamespace()/torch_available() branches for the torch/ggplot2/GGally/ggdensity/cli "package not installed" error messages (these only ever ran in an environment that actually lacked the package, which the coverage CI job never is), zero-row and zero-draw edge cases in the simulator pipeline, mcmc_init()'s two failure-to-start branches and split_rhat()/bulk_ess()'s degenerate-run guards, and npe_sequential()'s proposal-batch-cap warning and x_obs validation branches. .github/workflows/test-coverage.yaml now sets NOT_CRAN: "true", since a few tests (the future-multisession path in R/parallel.R, the slower two-moons NPE fit) skip_on_cran() to stay off CRAN's check machine, a concern that does not apply to this internal coverage job (#109) (#135).builtin_bar() (R/progress.R) previously extrapolated from elapsed / pos since the bar started, which is the average over the whole run so far. For NPE/NLE training that average is dragged down by the first epoch's one-time setup cost (device/net initialization, first CUDA/MPS kernel compilation, R JIT warmup) and, across n_restarts > 1, blends in restarts whose per-epoch cost can differ from the current one. A new bar_rate() estimates steps/sec from the last 8 calls instead, falling back to the lifetime average only while that window is too short to trust; bar_eta() turns the resulting rate into seconds remaining. The moving-target denominator that train_progress_total() projects (best_epoch + patience, growing every time validation loss improves) is unchanged and remains documented in ?nsbi_progress as a running estimate, not a promise -- this fix addresses the rate half of the ETA, not the total (#108) (#134).npe()/nle() gain a device argument for training on a GPU. device = "cpu" (the default, matching Python sbi; GPU is opt-in, never auto-selected), "cuda", "mps", or "gpu"/"auto" to resolve CUDA -> MPS -> CPU. Wrapping a fit in torch::with_device() used to crash with a mixed-device error: torch_tensor() built from an R matrix ignores that context and always lands on CPU, and MAF, MDN and NSF all build scratch tensors (torch_zeros(), an autoregressive rev_idx, the spline's tail padding, the MDN's Cholesky buffer) the same unqualified way, so a net moved to a GPU still hit CPU tensors partway through a forward or inverse pass. train_conditional_de() (R/train.R) now builds training/validation tensors on the resolved device and calls net$to(device = ); the estimator-specific forward/inverse code (R/flows.R, R/mdn.R, R/nsf.R) derives every scratch tensor from an existing tensor's device instead of relying on the default; and the de_log_prob()/de_sample() S3 methods, plus the MDN's separate i.i.d. fast path used by nle()'s MCMC sampling (R/likelihood.R) and stan_code()'s weight export (R/stan.R), move inputs to the estimator's device and results back to the CPU for the caller. "cuda"/"mps" error clearly when the requested backend is not available on this machine, rather than silently downgrading to CPU -- a request for a specific device that is not honored should not look like a slow CPU run; only "gpu"/"auto" falls back to CPU without complaint, since it never named a specific device. device is a no-op, not an error, for density_estimator = "linear_gaussian", which has no GPU concept. A fitted estimator remembers the device it trained on (fit$device, fit$de$device, stored as a plain string so it survives save_npe()) and load_npe()/load_nle() always rebuild the network on the CPU on reload, regardless of the device it trained on, since torch::torch_load() already defaults there. Closes #82 (#132).stan_fn_lingauss(), stan_fn_mdn() and stan_fn_maf() (R/stan.R) each wrote their own copy of the generated _sum_lpdf entry point -- the theta standardization, the x center/scale declarations, real total = 0, the for (n in 1:rows(x)) loop, and the jacobian-corrected return -- differing only in the loop body and in what, if anything, gets precomputed once outside the loop (the linear-Gaussian mean/Cholesky factor and the MDN's mixture head; MAF has nothing to hoist). R/stan.R now carries one internal helper, stan_sum_lines(fit, P, body, precompute = ""), that all three call instead of repeating the block; each generator keeps only its own loop body and precompute expression inline. Generated Stan output is unchanged for all three estimators. stan_sum_lines() sums the same per-observation log density that de_log_lik_iid() (R/likelihood.R) sums on the R side, so a comment now points from one to the other (#64) (#126).plot_coverage() and plot_tarp() each computed their own 99% Monte-Carlo binomial band by hand and built the same geom_ribbon() + geom_abline() + coord_equal() + theme_minimal() calibration figure around it, and plot_sbc() computed a third version of that band, unscaled, for its histogram's reference lines. R/plotting.R now carries two internal helpers, binom_band(nominal, n, level = 0.99) and calibration_plot(df, band, xlab, ylab), so the confidence level and the binomial parameterization live in one place instead of three: plot_coverage() and plot_tarp() call both and add only their own curve and title on top, and plot_sbc() calls binom_band() and scales it back to a count. All three functions keep their exact signatures and output (#63) (#125).resolve_x() (R/posterior.R) and resolve_x_iid() (R/nle_posterior.R) restated the same six lines -- fall back to the posterior's x_obs, stop if neither is supplied, check_numeric(), check_finite(), as_theta_matrix() -- differing only in whether row 1 is kept (with a warning) or every row is kept. R/posterior.R now carries one internal helper, resolve_obs(post, x, first_row, arg), that both call instead of repeating the body; resolve_x() and resolve_x_iid() keep their exact signatures and behavior, so no call site changed (#62) (#124).sbc()'s rank denominator now comes from the draws a trial actually returned, not from the requested n_posterior_samples. The chi-square uniformity test and expected_coverage() both bin sbc()'s ranks against a fixed scale, and that scale was read straight from the n_posterior_samples argument rather than from nrow(draws). The two numbers happen to always agree today: a short draw (a bounded prior and a leaky estimator defeating rejection sampling) is already turned into a hard error by diagnostic_draws() (0.4.3) before any rank is computed, so no trial reaches the binning step with fewer draws than requested. But the binning code was trusting an argument instead of reading back what actually happened, which is the wrong thing to trust even when the two currently coincide. The shared trial loop introduced below now reports the draw count it actually saw, and both sbc() and tarp() (whose own per-trial mean() was already dividing by nrow(draws), and so was unaffected) size their scale from that instead (#122).sbc() and tarp() each restated the same preamble -- check the fit and prior, draw truths from the prior, simulate, drop failed simulations -- and the same progress-bar/tryCatch trial loop, differing only in the metric computed per trial. R/diagnostics.R now carries two internal helpers: sbc_draws() for the preamble, which is also where the prior-width check now lives so both diagnostics get it from one place, and for_each_trial() for the loop, which both sbc() and tarp() call with their own f(draws, i) (a rank row for sbc(), a scalar coverage value for tarp()) instead of repeating either block (#61) (#122).Internal cleanup, no user-visible behavior change for print.nsbi_npe() and print.nsbi_nle(): both printed the same five blocks by hand -- parameter names, outcome names, the simulation count and any drops, the best validation loss, and the dead-network warning -- differing only in the save-function name (save_npe() vs save_nle()) and a "per observation" suffix on the data line. R/utils.R now carries two internal helpers, cat_fit_common() and cat_dropped(), that both print methods call instead of repeating the blocks; print.nsbi_sbc() and print.nsbi_tarp() also switch to cat_dropped() for their identical "N further trials dropped" line. print.nsbi_snpe() now calls cat_fit_common() too, which means a sequential fit's summary gains the parameter/data dimensions, dropped-simulation count, best validation loss and dead-network warning it never showed before, alongside the round table and not-amortized warning it already had (#60) (#121).
Internal cleanup, no user-visible behavior change: the MADE trunk (made_module() in R/flows.R, nsf_made_module() in R/nsf.R) and the plain MLP trunk (the MDN in R/mdn.R, the embedding network in R/embedding.R) each rebuilt the same linear, relu stacking loop. R/flows.R now carries one internal helper, mlp_layers(dims, masks = NULL), that both variants call -- masked when made_masks()$hidden is passed, plain otherwise -- so the four call sites share one loop instead of four (#59) (#120). The autoregressive forward/inverse transform-stack loops in flows.R/nsf.R are left untouched; the issue flagged those as a separate, riskier dedup.
fit_density_estimator() restated every fit_*() training default -- n_components, hidden, max_epochs, batch_size, lr, validation_fraction, patience, n_restarts, clip_grad_norm -- as a %||% fallback, once per estimator branch. Both callers, npe() and nle(), already pass every argument explicitly, so none of those fallbacks had ever fired; the numbers just had to be kept in sync by hand across fit_density_estimator(), each fit_*() signature, and the two callers. It now forwards only the arguments the chosen fit_*() accepts (intersect(names(dots), names(formals(fn))), then do.call()), so each fit_*()'s own signature is the only place its defaults live. linear_gaussian only takes ridge and verbose, so an n_components or embedding_net passed alongside it is still silently dropped, as before (#58) (#118).de_log_prob.* and two behind de_sample.* (de_sample.nsbi_de_mdn samples its mixture directly and shares nothing here), plus near-identical fit_mdn()/fit_maf()/fit_nsf() wrappers around train_conditional_de(). R/density_estimator.R now carries three internal helpers next to the generics they implement -- de_log_prob_torch(), de_sample_flow() and fit_torch_de() -- and each estimator's file calls them instead of repeating the body. fit_mdn(), fit_maf() and fit_nsf() keep their exact signatures; only their bodies collapsed, so man/fit_mdn.Rd, man/fit_maf.Rd and man/fit_nsf.Rd needed no edit (#57) (#117).slice_sample_run() has always counted n_evals, but nothing past slice_sample() itself read it. The count is attached as an n_evals attribute on the diagnostics data frame next to rhat and ess_bulk -- an attribute rather than a column, since it is one number for the whole run, not one per parameter -- and print() on an nsbi_nle_posterior now includes it alongside Rhat and ESS. It is the cost the adapted slice width is trying to keep down (see ?nsbi_mcmc), and there was previously no way to see it (#114).helper-torch.R's skip_if_no_torch() now calls the package's own torch_available() instead of re-implementing the same check inline, and require_torch() uses it too. drop_failed_sims() no longer returns an ok logical that nothing read. dmvnorm_chol() dropped its log argument -- every call site passed log = TRUE, so the exp() branch never ran. builtin_bar(), hint_parallel() and prior_scale() dropped default arguments no caller ever overrode. A doubled verbose guard in the training loop (R/train.R) is now a single cat() inside the if (verbose && ...) block that was already checking it. The commented-out neural-likelihood-estimation navbar entry in _pkgdown.yml is restored, now that the vignette it links to builds cleanly (#33) (#114)..github/workflows/release-check.yaml watches DESCRIPTION on main (plus a weekly schedule) and tags and publishes a prerelease once the version has had a minor or major bump, or once three weeks have passed since the last release tag on a patch-only bump; a patch bump inside that window publishes nothing (issue #110). .github/workflows/cran-submission.yaml fires on that prerelease and runs coatless-actions/cran-submission, which checks the package and submits the tarball. This release is the first tag the new workflow cuts, establishing the vX.Y.Z tagging habit CLAUDE.md has called for since 0.4.1 (#115).summary() now works on an nle() fit. Only summary.nsbi_npe was registered, so an NLE fit fell through to summary.default and printed a Length/Class/Mode table of the raw fit list, the de element holding the torch module included. Every other user-facing verb was extended to NLE in 0.4.3. summary.nsbi_nle calls the NPE method, which reads only fields both classes carry and dispatches print() on the object, so the NLE fit reports its data dimension per observation the way print() already does.plot_sbc() now takes a parameter name as well as an index, and checks whichever it is given. plot_sbc(res, param = 99) reached ranks[, 99] and reported subscript out of bounds, which names neither the argument nor how many parameters there are. The rank matrix carries colnames and every other plotting function labels by name, so param = "sigma" works too; a name that is not among the columns is refused with the ones that are.pairplot() checks limits against the number of parameters. A list or a matrix with the wrong number of entries indexed past its end and gave the same subscript out of bounds. The message now says how many limit pairs it got and how many parameters samples has.expected_coverage() requires levels strictly between 0 and 1. levels = c(-1, 2) used to be scored anyway: the central interval comes out empty or covers the whole line, so the row reads as coverage 0 or 1 and looks like a verdict on the fit. plot_coverage() inherits the check by passing levels through.print() on an NLE posterior reports unusable MCMC diagnostics as unavailable. split_rhat() and bulk_ess() return NA for a run with too few iterations, one chain, or a coordinate that never moved, and max(na.rm = TRUE) over all-NA is -Inf with a warning, so the line read max Rhat -Inf, min bulk ESS Inf. A partially scored run reports the parameters that did get a number and counts the ones that did not.stan_data() checks that the fit still has a live network, which stan_code() has done since it was written. A fit restored with plain readRDS() used to get the "save it with save_npe()" message from one and a dangling external pointer from the other, for the same fit and the same cause.load_npe() and write_stan_model() check the path they are given. load_npe() handed anything straight to readRDS(), which reports "cannot open the connection" and names the file only in the warning beside it; it now says the file does not exist, and refuses a value that is not one path the way save_npe() always has. write_stan_model() checks file before it transpiles the network, so a wrong destination costs no code generation.simulate_for_sbi() now recognises a call with simulator and prior the wrong way round. It takes the simulator first and npe(), nle() and npe_sequential() take it second, so the two functions a user calls in the same breath disagree. Getting it backwards used to fail inside sample_prior(), whose stopifnot() reports inherits(prior, "nsbi_prior") is not TRUE about a local that holds the simulator. A user reading that concludes the prior object is broken. inherits(simulator, "nsbi_prior") && is.function(prior) can only be the swap, since a prior is never a function and a simulator is never an nsbi_prior, so that case now says so and shows the right call. ?simulate_for_sbi says the order is reversed.simulate_for_sbi() also checks that simulator is a function and that prior is an nsbi_prior before it draws anything. npe_sequential() has had the simulator check since it was written; simulate_for_sbi() left both to surface further down.sbc() and tarp() now check that the prior they are given covers the parameters the fit was trained on. Both take prior = fit$prior as a default, so the argument is there to be overridden, and passing something else is legitimate: npe_sequential() trains on truncated proposals, and SBC against a narrower prior is a fair local calibration check. Nothing checked the width, though. sbc() sizes its rank matrix from fit$dim_theta and draws the truths from prior, so a mismatch reached sweep(), which recycles one against the other and ranks the truth against a comparison nobody asked for; in tarp() the same mismatch reached the z-scoring and the distances. Both now call check_prior(prior, dim = fit$dim_theta) before any simulation runs, so a wrong prior costs no simulator calls. ?sbc and ?tarp say what overriding prior changes about the diagnostic.c2st() now checks that its two sample sets are the same width, and that it has draws enough for the folds it was asked for. Two sets of different widths met each other at rbind(), which reports "numbers of columns of arguments do not match" and names neither x nor y; the message now gives both widths. More folds than draws was worse, because it was silent: rep_len() left the last folds empty, mean() of an empty test fold is NaN, and the accuracy came back NaN with no error. n_folds must now be fewer than the number of draws in the smaller set, and the message says how many draws that is. x and y also go through check_numeric(), so a character column is named rather than coerced to NA.c2st() assigned its cross-validation folds with the package's own sample() generic, which happened to work only because sample.default() forwards to base::sample(). It calls base::sample() now, as the subsampling ten lines above already did.posterior() took x_obs on trust, and an NA in it travelled through standardization into the density estimator and came back as all-NaN draws. The first complaint was stats::quantile()'s "missing values and NaN's not allowed if 'na.rm' is FALSE", raised from inside summary() and saying nothing about the observation. On the NLE side the same NA made every MCMC starting point non-finite, and the run failed reporting an initialization problem. The observation is the one input a user types out rather than simulates, so it is the one most likely to carry an NA from a real data set, and it was the only one the package did not check: drop_failed_sims() has always discarded non-finite simulations with a count and a rate, and slice_sample_run() has always refused a non-finite starting point. posterior() on either fit type, and an observation passed straight to sample(), log_prob() or map_estimate(), now go through check_finite(), which names the argument, counts the bad entries, says whether they are NA, NaN or Inf, and gives the position of the first one. It errors rather than warns: there is nothing sensible to condition on.resolve_x_iid() takes an arg argument so the message names the argument the caller used, obs for sample() and x for log_prob().as_sim_draw() has always reported everything wrong with a simulator's return value, down to which column was not numeric, but a failure that happens before there is a return value got none of that. npe(prior, function(mu, ls) rnorm(1, mu, exp(ls))) on an unnamed prior answered argument "ls" is missing, with no default, which is R's description of the symptom and says nothing about the cause: an unnamed prior sends the whole parameter vector to the first formal, so the second one never gets a value. Under a future plan the same error crossed a worker boundary before anyone saw it. The message is now Simulation 3 failed: <the original message>, followed by the parameter values that produced it and a pointer to ?nsbi_simulator. At most six values are shown, so a 40-parameter model does not fill the console. The original condition is kept as the parent and its class is carried onto the re-raised error, so a simulator that signals a custom condition can still be caught by class.call_sim_once() in R/simulator.R wraps the per-draw call, and describe_params() joins describe_value() in R/check.R.linear_gaussian estimator now checks the width of x like every other estimator. fit_linear_gaussian() recorded dim_theta and not dim_x, so its two methods had nothing to compare an incoming x against and passed it straight to the matrix product. An x of the wrong width came back as R's "non-conformable arguments", which names neither the argument nor the width expected of it, while the MDN, MAF and NSF all answer "Expected 3 columns but got 4". That gap mattered more here than it would anywhere else: linear_gaussian is the default in the examples, the torch-free path, and the oracle the test suite is built on. A bare vector of length dim_x is now also read as one observation rather than a column of values, again matching the neural estimators. Fits serialized before dim_x existed keep working.R/density_estimator.R gains a test file. fit_linear_gaussian(), lingauss_mean() and the two nsbi_de_lingauss methods were named nowhere in the suite despite being what everything else is checked against.prior_custom() now checks its arguments, and probes sample_fn and log_prob_fn once at construction. It is the one prior a user writes by hand and it took everything on trust, so the mistakes surfaced far from their cause. A log_prob_fn returning a single number instead of one density per row of theta passed the probe in nle_potential() and came back from the sampler as "Some MCMC starting points have zero posterior density. This is an initialization failure, not a sampling one", which is accurate about where it noticed and wrong about the cause. A lower or upper of the wrong length was recycled by sweep() inside within_support() with nothing but R's "STATS is longer than the extent of 'dim(x)[MARGIN]'" warning, and the wrong support test then decided which posterior draws were rejected as leakage and what log_prob() renormalized by. A sample_fn whose width disagreed with dim was caught, but anonymously, as "Expected 2 columns but got 1". dim must now be a positive whole number, sample_fn and log_prob_fn functions of one argument, and lower/upper numeric of length dim with upper above lower. sample_fn(2) is called once at construction and has to return a 2 x dim numeric matrix, and log_prob_fn is evaluated on those two rows and has to return two numbers. Every message names the argument it rejected.lower and upper accept a single number, recycled to every parameter, so a positive-support prior is lower = 0 rather than lower = rep(0, dim). The recycling is deliberate and documented; it is only the silent recycling of a wrong-length bound that is gone.prior_custom() gains param_names. new_prior() has always carried parameter names, and prior_uniform()/prior_normal() take them from the names of low/mean, but prior_custom() had no way to pass them. That was not cosmetic: sim_dispatch() decides from prior$param_names whether a simulator is called with one scalar per formal or with the whole parameter vector, so a custom prior could never use the named-simulator signature. ?prior_custom also now says why a custom prior cannot be written out by stan_code().print() on a prior bounded on one side only no longer errors. It printed upper whenever lower was set, and prior_custom(..., lower = 0) with no upper is a normal thing to build.check_function() and check_bound() join the validators in R/check.R.npe() and nle() now warn when a column of theta or x has no spread. Standardization has always guarded against dividing by zero by leaving such a column at scale 1, and it did so without saying anything. A constant summary statistic does not respond to the parameters, so it tells the estimator nothing and the mistake is invisible from the outside: training converges and the posterior looks plausible while the coordinate does nothing. The warning names the columns, by name where the matrix has column names and by index otherwise, and says which side of the fit they are on. A single row gets its own wording, since sd() of one value is NA rather than zero. The standardize = FALSE path stays silent: its degenerate standardizer is built from a one-row zero matrix on purpose.fit_standardizer() takes a what argument naming the side it is standardizing, the way drop_failed_sims() already did. Callers that leave it NULL warn about nothing.theta or x is now an error that names the column. The pre-computed path coerced its input with storage.mode(x) <- "double", so a character or factor column turned into a column of NA with R's generic "NAs introduced by coercion" warning. Every row was then dropped as non-finite and the run stopped with "All 50 simulations returned non-finite output (NA, NaN or Inf). Check the simulator on a single prior draw", which is wrong twice over: no simulator ran, and the data was fine apart from one text column. npe(prior, theta =, x =), nle() and log_prob() now say `x` has non-numeric columns: b, the way as_sim_draw() has always reported the same mistake in simulator output.drop_failed_sims() no longer tells you to check the simulator when you did not call one. On the pre-computed path it points at theta or x instead.check_numeric() carries the type half of check_matrix(), so entry points that disagree about shape (a bare vector is one parameter set to log_lik() and a column of values to npe(theta =, x =)) still share one message about types.n_simulations = -5 failed inside stats::runif() as "invalid arguments", validation_fraction = 1 inside seq() as "wrong sign in 'by' argument", batch_size = 0 as "invalid '(to - from)/by'", and n_restarts = 0 skipped the restart loop and came back as "Training failed: no restart produced a finite validation loss", which blames training for an argument. n_simulations = 0 reached cbind() and warned about recycling; n_simulations = 1 trained on one draw and said nothing. The counts (n_simulations, n, n_sbc, n_tarp, n_posterior_samples, n_folds, n_rounds, n_init, n_normalization, n_truncation_samples, max_sampling_batches, max_proposal_batches), the training controls (max_epochs, batch_size, lr, validation_fraction, patience, n_restarts, clip_grad_norm) and the architecture arguments (n_components, n_transforms, hidden, n_bins, tail_bound) are now checked at the top of the entry point that takes them, in npe(), nle(), simulate_for_sbi(), npe_sequential(), sample(), log_prob(), map_estimate(), sbc(), tarp() and c2st(). Each message names the argument and the value it rejected. npe() and nle() resolve density_estimator there too, so a misspelled estimator is an error before the simulator runs rather than after.validation_fraction is also checked against the number of simulations, inside train_conditional_de(), which is the first point where that number is known. A fraction that leaves no training rows now says how many rows the split would need instead of failing in seq().check_counts() validates a vector of counts (hidden), and check_positive() gains allow_inf for clip_grad_norm, where Inf is the documented way to disable clipping.R/check.R: check_matrix(), check_count(), check_prob(), check_positive(), check_prior() and check_finite(). Each names the argument it rejects and says what was actually wrong, following the voice of the simulator-output checks. They are internal, and the entry points adopt them one at a time.log_lik() now rejects a wrong-length theta or x instead of reshaping it. A vector shorter or longer than the fit's width was recycled into a matrix of the right number of columns, so log_lik(fit, c(0.1, 0.2, 0.3), x) on a two-parameter fit returned two log-densities computed from recycled values, and a length-4 vector became two parameter sets with no message at all. At a public boundary a bare vector now means a single row, or it is an error. A matrix whose row count matches the expected width is told to transpose. as_theta_matrix() keeps its permissive reshaping for internal callers that rely on it.resolve_x() kept row 1 and said nothing about the rest, which is right for NPE and invisible to the user. The trap is the asymmetry with nle(), where the rows of x_obs are independent observations that the log-likelihood sums over: the same matrix means "200 observations" to nle() and "the first observation" to npe(), so moving a working call from one to the other silently dropped all but one data point. The warning names the row count and points at nle(). It stays a warning because taking row 1 of a simulation matrix is a reasonable thing to ask for, and sbc()/tarp() pass single rows anyway.nle(). Where npe() learns the posterior directly, nle() learns a surrogate likelihood \eqn{q(x \mid \theta)}. The reason to want that is repeated observations. An NPE fit is trained for one fixed data dimension, so conditioning on \eqn{n} independent trials means retraining for every \eqn{n} or compressing to summary statistics; NLE learns the density of a single trial, so the log-likelihood of \eqn{n} of them is a sum and \eqn{n} is free at inference time. The cost is that the posterior is no longer a forward pass. For one fixed observation with high-dimensional data, npe() is still the better choice, and ?nle says so.log_lik(fit, theta, x) evaluates the surrogate likelihood, summing over the rows of x as independent observations, and likelihood_fn(fit, x_obs) returns it as a plain vectorized function(theta). That closure is the point of contact with the rest of R: it goes straight into optim(), an MCMC package, an importance sampler, or a profile likelihood, with nothing downstream needing to know about neuralsbi.posterior() on an nle() fit returns an MCMC-backed posterior and gains sampler, n_chains, warmup, thin and init_strategy. The default sampler is a vectorized univariate slice sampler, matching Python sbi's default: nothing to tune, no dependency, and bounded priors need no special handling. Its vectorization runs across chains, so one step costs one batched forward pass rather than n_chains separate ones. The slice width adapts to the target during warmup, which matters because a width inherited from the prior is badly wrong once many observations have concentrated the posterior, and every unit of mismatch is paid for in wasted density evaluations. Draws carry split-Rhat and bulk ESS, and are cached on the posterior so summary() and repeat sample() calls do not re-run a chain. n_chains defaults to 20 under "slice" and 4 under "stan", because a slice chain rides along in a batch someone else is already paying for and a Stan chain is a process with its own warmup.thin defaults to 2. With the adapted width, thin = 2 already gives a bulk ESS of about 96% of the retained draws on a Gaussian target, so thinning harder buys very little for what it costs, and the reported ESS is there to tell you when to raise it. Python sbi thins by 1 (it thinned by 10 up to v0.21), so this sits between the two.theta to a Gaussian mixture over x and never sees x, so for n observations the network runs once and all n densities come off the same mixture; the linear-Gaussian baseline behaves the same way. A flow's transforms depend on x too, so it has to run n times. With a few thousand observations that is the difference between seconds and minutes per MCMC step, and it is worth weighing when choosing an estimator for repeated data.n_theta x n_obs matrix, the largest object in the loop, is never built. And an MDN's chunking is sized by the pair count rather than by pairs times components, which had been splitting a 5000-observation call ten ways for a 400 KB intermediate. On the four-parameter g-and-k model in vignette("neural-likelihood"), 2000 draws from 20 chains went from 117s to 25s at 500 observations and from 342s to 68s at 5000, with split-Rhat and bulk ESS unchanged.torch::jit_trace() records the same code and replays it in one crossing. It is a shortcut, not a path. Nothing is recorded until an evaluator has been called a few times, so a single log_lik() never pays for a compiler; a trace fixes the shapes it saw, so there is one per parameter-row count and each is checked against the eager result before anything uses it; and tracing or checking failing just leaves the eager path in place. options(neuralsbi.jit = FALSE) turns it off.log_prob() on an NLE posterior returns the unnormalized log posterior. The evidence is not available, so normalize is ignored with a warning rather than returning a number that looks normalized and is not.stan_code(), stan_data() and write_stan_model() export a fitted likelihood as Stan source. The generated functions block recomputes \eqn{\log q(x \mid \theta)} in Stan's own language with the trained weights passed as data, so Stan differentiates it and NUTS gets exact gradients, with nothing linked against torch at run time. The result is an ordinary Stan function of theta, which is what makes it worth having: the surrogate stops being the whole model and becomes one term in a model you write, next to a hierarchical prior, covariates, or a second data source whose likelihood you do know. "mdn", "maf" and "linear_gaussian" are supported; "nsf" is refused with a message naming the alternatives. posterior(fit, x_obs, sampler = "stan") runs the generated model through cmdstanr or rstan and returns draws like any other path.sbc(), tarp() and posterior_predictive() accept an nle() fit, and sbc()/tarp() forward ... to posterior() so the MCMC controls reach it. Every SBC trial is a separate MCMC run, so start small.save_npe()/load_npe() handle nle() fits too, with save_nle()/load_nle() as aliases.npe() and nle() now share prepare_simulations() instead of each carrying its own copy of the simulate/coerce/drop/standardize preamble.sample_posterior() now goes through the sample() generic instead of calling sample.nsbi_posterior() directly. An nle() posterior inherits nsbi_posterior, so the old call ran the NPE forward-pass sampler against an estimator whose target and condition are swapped. That errored with non-conformable arguments when dim_theta and dim_x differ, and returned draws from the wrong distribution without complaint when they happen to match.posterior() on an nle() fit now validates n_chains, warmup and thin instead of coercing them with as.integer(). thin = 0 ran the warmup and kept nothing, so the slice sampler returned its zero-initialized array and every draw was 0, with no error and an Rhat computed on those zeros. Values below the bound, fractional values and NA are now errors.theta as well as x. A row was dropped only when its simulator output was non-finite, so an NA among pre-computed parameters passed straight through to the estimator: npe(prior, theta = theta, x = x) with one missing parameter failed inside chol() with "the leading minor of order 1 is not positive", and on the torch path it became "Training failed: no restart produced a finite validation loss", which blames training for a bad input. Such rows are now dropped like any other failure, and the warning says which side was non-finite ("parameters", "output", or "parameters or output").npe_sequential() now checks x_obs against the simulator's output width, at the end of round 1, which is the first moment that width is known. An x_obs of the wrong length was reshaped into several rows and only the first was targeted, so the run truncated its proposals around a value the caller never gave while print() reported all of them as targeted. A TSNPE fit is not amortized, so targeting the wrong observation is the whole failure. A missing, NULL, non-numeric, NA or multi-row x_obs is an error as well, and n_rounds must be a single integer of at least 1: n_rounds = 0 skipped the loop and returned a bare list classed nsbi_snpe.sbc() and tarp() now error when a trial's posterior returns fewer draws than n_posterior_samples instead of scoring it anyway. sample() comes back short when a bounded prior and an estimator that leaks mass outside it defeat rejection sampling, and it only warns. sbc() then binned those ranks against n_posterior_samples while they were drawn from a smaller set, which compressed every rank toward zero: the chi-square uniformity test rejected and expected_coverage() fell below the diagonal, so lost draws were reported as a miscalibrated posterior. Rescaling the short trial on its own is no better, since it is then scored on a different resolution from the rest, so both functions stop and name the trial and the shortfall.npe() and nle() now match density_estimator against the allowed names before the simulator runs. The check lived in fit_density_estimator(), which is reached only after the simulations are in hand, so density_estimator = "mfa" spent the whole budget and then failed on the typo. For the expensive simulators this package is for, that is minutes to hours. The same hoist fixes the embedding warning, which compared the unmatched value: match.arg() accepts abbreviations, so density_estimator = "linear" selected linear_gaussian and dropped embedding_net without saying so. fit_density_estimator() keeps its own match.arg() because it is also reachable directly.prior_uniform(c(a = -3, b = -1), c(a = 3, b = 1)) with function(mu, ls = 0) sent the whole two-parameter vector to mu, left ls at its default, and trained on the result without complaint. One typo in a prior name was enough. The warning names the parameters that found no formal and points at ?nsbi_simulator. It is a warning and not an error because a vector-signature simulator whose first argument happens to carry a parameter's name is legal.pairplot()'s lower triangle now shows highest-density regions (via the new ggdensity Suggests, geom_hdr(), at the 50/80/95/99% probability levels) instead of a raw scatter of draws. A scatter of 10,000+ points overplots into an undifferentiated blob at the resolution most posteriors are viewed at; nested HDR contours show where the mass actually concentrates. The truth cross-hair markers are unchanged. col now sets the region fill colour; alpha applies only to the diagonal marginal densities, since the lower triangle shades itself by probability level.pomp or deSolve -- so meeting the old contract meant wrapping the model in an apply loop and thinking about column-major indexing before anything ran. Two signatures are accepted, decided from formals(simulator): when every parameter name appears among the formals the parameters arrive by name, one scalar each (function(mu, sigma) ...); otherwise the whole named parameter vector goes to the first argument (function(theta) ...). A simulator returns a numeric vector of length d, a scalar, or a one-row matrix or data frame, and names on that output become the outcome names. See ?nsbi_simulator for the contract and the migration examples.simulate_for_sbi(function(theta) c(theta[1] * 2, theta[2] * 2), prior, 100) returned a 100 x 1 matrix of nonsense with no error and no warning; it now returns the 100 x 2 matrix it should.npe(), simulate_for_sbi(), npe_sequential(), sbc(), tarp() and posterior_predictive() gain sim_args, a named list forwarded to every simulator call. Observed data, a time grid, a fixed population size, a design matrix or a solver tolerance travel from the call site instead of being captured in a closure -- which also keeps them out of what gets serialized to a future worker. A list rather than ... because every one of npe()'s formals sits before ..., so R's partial matching would silently capture x, theta, n or seed.NA, NaN or an infinite value are dropped, together with their parameters, with one warning per run reporting the count and the rate. A single non-finite value would otherwise poison the training loss and surface much later as a NaN validation loss. The count is recorded on the fit and shown by print(). Nothing left is an error. In sbc() and tarp() a failed draw removes the whole trial; in posterior_predictive() it reduces the number of predictive draws. Pre-computed theta/x passed to npe() are checked the same way. Note that dropping conditions on the simulator having succeeded: when failure depends on the parameters, the fit targets the posterior given success, so a high drop rate is a modelling signal.future plan and whatever the worker count. The 0.4.0 guarantee held only at a fixed chunk size.chunk_size is gone, from npe(), simulate_for_sbi(), npe_sequential(), sbc(), tarp() and posterior_predictive(), along with options(neuralsbi.chunks). Chunking existed in 0.4.0 to make results reproducible across backends: the split had to depend on n alone, which meant users had to know about it and could change their draws by changing it. Per-simulation RNG streams give that guarantee outright, so what is left is ordinary parallel scheduling. Batches are now sized from the worker count, cannot affect a result, and are not a setting. Running sequentially there are no batches at all, just a loop.n_bins and tail_bound are explicit npe() arguments; npe() no longer takes ....save_npe() / load_npe(). A fit whose estimator is "maf", "mdn" or "nsf" holds a torch module, which is an external pointer: saveRDS() writes the pointer, the file reloads without complaint, and the first call that touches the network fails with external pointer is not valid. save_npe() writes the weights with torch::torch_save() and everything else as ordinary R objects, into one .rds; load_npe() rebuilds the network from the recorded architecture and restores them. An overnight fit can now be reloaded the next morning, which is what amortization was for. posterior() and print() also detect a fit that came back from readRDS() and say so, instead of failing later with a torch error.future plan -- library(future); plan(multisession) -- and every function that calls a simulator (npe(), simulate_for_sbi(), npe_sequential(), sbc(), tarp(), posterior_predictive()) spreads the work across workers. There is no new argument to pass and no parallel variant to call: with no plan declared everything runs sequentially as before, and neuralsbi mentions the two lines above once per session (options(neuralsbi.parallel_hint = FALSE) to silence it). Each chunk of parameters draws from its own L'Ecuyer-CMRG stream, so a given set.seed() produces the same simulations sequentially and on any number of workers. See ?nsbi_parallel.progressr installed, neuralsbi emits standard progressr updates, so progressr::handlers() and with_progress() control reporting; without it, a built-in bar needing no extra packages does the job. The training bar targets the epoch at which early stopping would fire and revises that target as the validation loss improves. See ?nsbi_progress.npe(), simulate_for_sbi(), npe_sequential(), sbc(), tarp(), and posterior_predictive() gain a chunk_size argument controlling how many parameter rows go to the simulator per call. The default splits a run into about 64 chunks; because the split depends only on the number of simulations, results do not change with the number of workers. A simulator whose output must be produced in one call can set chunk_size to the full simulation budget.future and progressr are Suggests, not dependencies; parallel (base R) moves into Imports.ggplot2/GGally are not installed, instead of failing. Both are Suggests, so R CMD check under _R_CHECK_FORCE_SUGGESTS_=false -- the configuration CRAN uses on a machine without them -- previously hit 5 errors from require_ggplot2(). The new skip_if_no_ggplot2()/skip_if_no_ggally() helpers mirror the skip_if_no_torch() contract already used for the neural tests, so the suite runs everywhere.vignette("sir-time-varying-beta") reworks the epidemic model so that its posterior predictive tracks the observed case peaks instead of overshooting them. The introduction day and seed size are now inferred per state rather than fixed at two infections on 2020-01-21 (which left the epidemic's phase to demographic noise: replicate simulations at one fixed parameter differed by a factor of several thousand at the peak bin), the ascertainment prior is pinned by spring-2020 seroprevalence (case counts alone identify only the product of ascertainment and incidence, and the unconstrained fit slides toward vast, barely-ascertained epidemics), and the regime-duration parameterization drops a coordinate that the simulator's rescaling had left unidentified. The article reports the effective reproduction number R~e~(t) = beta(t) S(t) / (N gamma) computed from the simulator's own susceptible trajectories, adds a prior-predictive check before training, and summarizes the posterior predictive by its median rather than its mean.prior_uniform()'s low/high or prior_normal()'s mean (e.g. c(beta = 0, gamma = 0)), or attach colnames() to a simulator's output, and those names carry through npe(), sample(), map_estimate(), sbc(), expected_coverage(), and posterior_predictive() without any extra arguments. plot_sbc(), plot_coverage(), pairplot(), and plot_posterior_predictive() use the names for titles, legends, and facet labels; a name that happens to be valid R syntax ("beta[1]", "rho", "sigma^2") renders as its plotmath symbol (Greek letters, sub/superscripts) instead of literal text.pomp package, vignette("sir-epidemic"). Both methods fit the same stochastic SIR epidemic: pomp via particle-filter MCMC (pmcmc), neuralsbi via neural posterior estimation. The vignette contrasts what each needs from the model — pomp a measurement density, neuralsbi only a simulator — overlays the two posteriors, scores their agreement with a C2ST, and confirms the neural fit with SBC. The comparison is precomputed, so pomp is needed only to regenerate the article, not to build or check the package.ggplot2 and GGally::ggpairs() instead of base graphics. pairplot(), plot_sbc(), plot_coverage(), plot_tarp(), and plot_posterior_predictive() keep their signatures (pairplot() gains an alpha argument) but now build and print a ggplot/ggmatrix object, returned invisibly for further customization. ggplot2 and GGally move to Suggests, following the same graceful-degradation pattern as torch: call any plotting function without them installed and you get an informative error, not a crash.vignette("intro-to-sbi"): a short beginner tutorial covering the three ingredients (prior, simulator, observation), amortized training, and a first calibration check, using a g-and-k distribution simulator whose likelihood has no closed form.README.Rmd, so the usage example runs at render time and its output cannot drift from the code. The example is a plain linear regression: the posterior recovers the ground-truth coefficients, which the reader can check against ordinary least squares.summary() and as.data.frame()
methods into the single summaries help topic (they had drifted into
separate .Rd files with duplicated \alias entries, which also produced
duplicate HTML anchors). Wrapped theta_{<d} in \eqn{} in the
made_masks docs so the Rd no longer drops braces. Dropped the bare "NPE"
acronym from the DESCRIPTION to avoid the spurious misspelling note.Defaults now match Python sbi, so a workflow reads the same in both
packages and results can be cross-checked. Changes to npe() defaults:
the density estimator is now "maf" (was "mdn"); MDN mixture components
default to 10 (was 5); NSF spline bins default to 10 (was 8); the training
batch size is 200 (was 100). max_epochs is raised to 2000 as a guard cap
that early stopping (patience = 20) normally reaches first, mirroring
sbi's effectively-unbounded epoch budget. lr, validation_fraction,
patience, clip_grad_norm, n_transforms, and hidden already matched.
Pass any of these explicitly to recover the previous behavior.
First CRAN submission. Dropped the development .9000 version suffix,
removed the redundant Author/Maintainer fields (now derived from
Authors@R), and tidied the package title.
Embedding networks (roadmap v0.4). embedding_mlp() builds a learned
summary network that maps raw observations to a low-dimensional feature
vector; pass it to npe(..., embedding_net = ) and the MDN, MAF, and NSF
estimators condition on the features instead of the raw data, training the
embedding jointly. The estimators still take raw x at the de_* boundary
(dim_x is unchanged), so sampling and log_prob route through the
embedding automatically. Ignored, with a warning, by linear_gaussian.
vignettes/<name>.Rmd.orig, and
vignettes/precompute.R bakes it into a static vignettes/<name>.Rmd
(results, printed values, and figures inlined). CI and pkgdown re-render
that static Markdown with no torch at build time, so the expensive neural
training runs once, locally, instead of on every build. Re-run
Rscript vignettes/precompute.R after editing any .Rmd.orig.inst/benchmarks/two_moons_calibration.R):
SBC, expected coverage, and TARP for a two-moons NSF fit, with figures
written to docs/figures/ (roadmap milestone M2).npe_sequential()). Removed a truncated duplicate
of the SIR vignette.DESCRIPTION (Pedro Nascimento de Lima, with ORCID).npe_sequential(): multi-round NPE targeting a single observation via
truncated-prior proposals (TSNPE, Deistler et al. 2022). Each round truncates
the prior to the highest-probability region of the current posterior and
retrains on all accumulated simulations; the standard NPE loss stays valid,
so no importance correction is needed. Returns an nsbi_snpe fit that works
with posterior(), sample(), and the diagnostics, but is only valid at the
targeted x_obs. Verified against the analytic linear-Gaussian posterior.tarp() diagnostic and plot_tarp() (Lemos et al. 2023): a joint
expected-coverage test using random reference points, complementing the
per-parameter sbc() ranks. Detects posteriors with calibrated marginals
but wrong correlation structure.plot_posterior_predictive(): marginal predictive histograms with the
observation marked; returns the observation's predictive quantiles.log_prob() integrates to one over the support and returns -Inf outside
it (test-posterior-normalization.R).R CMD check failed on three counts: the npe() example required
libtorch (it now uses the torch-free linear_gaussian estimator and runs
unconditionally), the hand-maintained npe.Rd/fit_mdn.Rd usage sections
had drifted behind the code (missing n_restarts, clip_grad_norm,
n_transforms, and the "maf"/"nsf" options), and CLAUDE.md was not in
.Rbuildignore. The test-torch job also failed because torch 0.17 refuses
a TORCH_HOME that does not exist; the workflow now creates it first.train_conditional_de()):
best-of-n restarts, learning-rate decay on plateau, gradient clipping,
per-epoch loss history.density_estimator = "maf") and Neural Spline
Flow ("nsf", autoregressive rational-quadratic splines) join the MDN and
the closed-form linear_gaussian baseline.task_gaussian_linear(), task_two_moons(),
task_slcp(), task_sir()) shared between tests and the
inst/benchmarks/ head-to-head benchmark harness.summary() methods, as.data.frame() tidy accessor, plot_coverage().R CMD check plus a test-torch job with cached libtorch.npe(), linear_gaussian
and MDN estimators, posterior sampling with leakage correction, SBC, expected
coverage, C2ST, posterior-predictive checks, pairplot(), plot_sbc().