Bootstrap/bagging confidence intervals (boot-sp merge PR-2: statdist + smoothed methods) - #818
Draft
DLWoodruff wants to merge 18 commits into
Draft
Bootstrap/bagging confidence intervals (boot-sp merge PR-2: statdist + smoothed methods)#818DLWoodruff wants to merge 18 commits into
DLWoodruff wants to merge 18 commits into
Conversation
Add mpisppy/confidence_intervals/bootsp/statdist/, a trimmed port of the statdist library holding only the univariate distributions the smoothed bootstrap methods need (uniform, normal, student-t, Gaussian kernel, epi-spline, empirical, discrete) plus their support modules. The multivariate machinery (copula.py, vine.py, bicop.py and the multivariate distribution classes) is dropped; that also removes the `from scipy.stats import mvn` import (gone in scipy 1.14) and the gosm hook. scipy is imported lazily via pyomo.common.dependencies so the empirical bootstrap path stays scipy-free, and matplotlib is imported inside the plotting methods so it stays optional. distribution_factory now scans the registry once. .ruff.toml relaxes the legacy-style rules this ported library predates, scoped to the statdist subpackage. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add smoothed_boot_sp.py (smoothed bootstrap and bagging that fit a statdist univariate distribution to the sampled data and resample from it) with compute_smoothed_ci as the single smoothed-dispatch point, the counterpart to boot_sp.compute_ci for the empirical methods. Wire the drivers to route smoothed vs empirical by boot_method: user_boot returns the empirical 6-tuple or the smoothed (ci_gap, center_gap) pair; simulate_boot gains smoothed_main_routine. This fixes the design-doc section 4.3 latent bug: boot-sp's smoothed simulation called the undefined fit_resample_utils.compute_xhat; it now calls boot_utils.compute_xhat. The PR-1 "smoothed not yet merged" error is removed. Non-root ranks return a matching-arity (None, None) so callers can unpack safely under mpiexec. The only scipy use (norm.ppf) is replaced by statistics.NormalDist().inv_cdf, matching boot_sp.py. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add the three examples that build their scenario data with statdist and so could not ship in the statdist-free PR-1: farmer (crop yields perturbed by a univariate distribution), cvar (Lam & Qian, standard normal data), and multi_knapsack (Vaagen & Wallace, deterministic data from a json file). Each is self-contained: a fixed-name xhat_generator builds the EF from the module's own scenario_creator (no external amalgamator coupling), with an empirical json/bash and a smoothed json. multi_knapsack's data_sampler reads its deterministic data from the file when the smoothed driver has not stashed it, so the empirical path works too, and resolves the json relative to the module directory. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add test_boot_sp_smoothed.py: direct tests of the statdist univariate distributions, the empirical farmer/cvar tests that PR-1 could not host (they need statdist), and the smoothed methods (kernel/bagging serial and under mpiexec -np 2; epi-spline tests skip without a nonlinear solver). Wire it into the confidence-intervals CI job and run_coverage.bash, serial and np=2, in the same commit. Update the two PR-1 test files whose "smoothed not yet merged" assertions no longer hold. Register a farmer (empirical) and a cvar (smoothed) run in part 1 of run_all.py. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the "smoothed methods merged separately" note with a description of the empirical vs smoothed families, add the five Smoothed_* tokens to the methods table, document data_sampler and the smoothed-only options, and add a section on the smoothed methods, the bundled statdist library, and the farmer/cvar/multi_knapsack examples. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The Test_empirical_examples tests inspected compute_ci's result on every rank, but only rank 0 gets the values (non-root ranks get a tuple of Nones). Subscripting None raised on the non-root rank and aborted the test mid-loop, so that rank skipped the next compute_ci while rank 0 entered its collective and deadlocked. Guard the value checks on rank 0 (and assert the None-tuple off-rank) so both ranks always make the same collective calls. Found by running the suite under mpiexec -np 2. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ning) The Gaussian-kernel distribution's pdf returned gaussian_kde.evaluate(x), a size-1 numpy array. The base-class cdf feeds pdf to scipy.integrate.quad, which converts the integrand to a scalar and, on NumPy>=1.25, emitted an "array to scalar" DeprecationWarning tens of thousands of times during a smoothed run. Return float(...[0]) so the integrand is a scalar; the numerics are unchanged. This removes the need for the warnings filter that was in test_boot_sp_smoothed.py. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Two defects in the smoothed bootstrap, both measured against Chen & Woodruff (2024), "Distributions and Bootstrap for Data-based Stochastic Programming". They are fixed together because the second is what forces the index-space reservation the first introduces. 1. smoothed_resample_helper advanced the record index by one per batch while taking a block of subsample_size consecutive records, so consecutive batches shared all but one of their draws -- a sliding window rather than independent resamples. Algorithm 3 draws a fresh set of N points from the fitted distribution for each of the B batches. Stride by the batch size, as smoothed_bagging already did, so the blocks are pairwise disjoint. The estimated spread was badly understated: on the cvar example (N=20, nB=10) the interval widens about fivefold. The dead boot_cfg copy goes with it -- it was mutated but never passed to the solves, which is what hid the missing stride. 2. smoothed_bootstrap called center_smoothed while use_fitted was still False, so the center was the purely empirical gap. Algorithm 3 takes the center from Algorithm 2 run on the fitted distribution, and that smoothed center is the paper's leading conclusion. Set use_fitted before estimating the center (smoothed_bagging already did), and reserve the center's block ahead of the batch blocks, since both now sample the same distribution and must not reuse each other's draws. simulate_boot's coverage replications are spaced to match. The ported spacing was a hard-coded nB * 100, unrelated to how many record numbers a replication actually consumes; it is now exactly that footprint. The empirical harness needs no spacing at all, because it gets independent numpy streams from default_rng([seed_offset, word]); the smoothed path addresses draws by record number, since that is what the model seeds each draw with, so its replications are separated by giving each a disjoint block. The new test captures the actual scenario pools and asserts the batch blocks are pairwise disjoint and disjoint from the center's, and that use_fitted holds when the center is computed; it fails on the old code for both reasons. No locked values change: the smoothed tests assert only structure. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
np.var defaults to ddof=0, but the algorithm specifies sample variances and the empirical estimators already use ddof=1 for their Gaussian half-width. ddof=0 understates s1 by (B_I-1)/B_I -- a third at the B_I=3 used in run_all.py -- so the bagging interval came out too narrow. Since s1 is the variance among the B_I per-seed-point averages, B_I < 2 leaves it undefined and was silently producing a nan; that is now a clear error. center_smoothed took an mpicomm it never used; dropped, along with the argument at its one call site. run_all.py's do_one_boot still described its examples as statdist-free with the others "merged separately"; it registers farmer and cvar now. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The estimators form every gap as (value at xhat) minus the optimal. For a maximization the value at xhat is a lower bound on the optimal, so that quantity is non-positive -- and user_boot floors the reported interval with ci_gap[0] = max(0, ci_gap[0]). A maximization run would therefore not merely be wrong, it would report [0, 0]. Per the repo-wide rule that maximization either works or raises, this is the raise, checked in solve_routine, which every extensive form goes through. The message says what supporting maximization would take rather than just refusing. The test builds one trivial model twice, once each sense, and checks that the maximization raises while the minimization gets past the guard and fails later on a deliberately bogus solver name, so the guard is not what stopped it. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Codecov Report❌ Patch coverage is Additional details and impacted files@@ Coverage Diff @@
## main #818 +/- ##
==========================================
+ Coverage 76.49% 76.57% +0.08%
==========================================
Files 175 183 +8
Lines 23236 24283 +1047
==========================================
+ Hits 17774 18595 +821
- Misses 5462 5688 +226 ☔ View full report in Codecov by Harness. 🚀 New features to boost your workflow:
|
Writing direct tests for the trimmed statdist univariate distributions turned up eight things that were wrong in the ported code: - Parameter.instantiated was inverted: it was True exactly when the parameter had no value. set_value now keeps it up to date too. - memoize_method left the keyword *values* out of the cache key, so cdf(x, epsabs=1e-4) and cdf(x, epsabs=1e-9) shared one answer. - memoize_method also died with "unhashable type: 'list'" when a method was handed a list or a dict, hiding the method's own behavior (region_expectation and region_probability validate their argument and say what is wrong). Such a call cannot be cached, so it is now passed through; the docstring, which claimed the arguments were converted, is corrected to match. - UnivariateDiscrete.var had the sign backwards (mean**2 - E[X**2]), so the variance of every non-degenerate discrete distribution came out negative. - UnivariateDiscrete accepted breakpoints whose probabilities summed to less than one: the check was one-sided. - UnivariateDiscrete.cdf_inverse's error message named cdf instead. - UnivariateGaussianKernelDistribution.fit dropped the bw_method it was given and always passed None, so the caller could not widen the bandwidth. - The kernel distribution's draw method passed [n, seed] to gaussian_kde.resample as the sample size, which raises; the size and the seed are two arguments. It is also named generates_X now, like its siblings. - UnivariateEmpiricalDistribution.cdf and cdf_inverse raised IndexError on a constant sample (a one-record sample included) while looking for a second point to extrapolate along. Every quantile of a constant sample is that value, and the cdf steps from 0 to 1 there. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
statdist arrived with the smoothed methods as its only exercise, which left it at 27% coverage. These are direct tests, all solver-free, of what the library promises: the plotting-position contract of the empirical cdf and its inverse, the closed forms for the uniform, normal and Student's t, the kernel density's padded domain and scalar pdf, the discrete distribution's moments and step function, and the registry's naming rules. The generic machinery in base_distribution.py (numeric cdf by quadrature, cdf inversion by refinement, expectations, sampling) is tested against a local distribution with a closed form for all of it -- density 2x on [0, 1] -- so the assertions are exact rather than self-referential. utilities.py gets tests for both memoizers and the argv context manager. Coverage of the reachable statdist code: distributions.py 37% -> 87%, utilities.py 48% -> 98%, base_distribution.py 24% -> 42% (only three univariate statements are left uncovered there: two abstract method bodies and the interactive plt.show(); the rest of the miss is the MultivariateDistribution class, which the trim left behind and which nothing can reach). Everything lands in the existing test file, so no harness wiring changes. Two tests are skipped in the usual environments: the plot smoke test needs matplotlib, and the check that importing statdist does not import scipy needs to spawn a plain python subprocess, which cannot be done from inside an mpiexec launch. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The port dropped copula.py, vine.py, bicop.py and the multivariate distribution classes in distributions.py, but MultivariateDistribution itself stayed behind in base_distribution.py along with the five decorators that exist only to serve it (fit_wrapper, accepts_dict, returns_dict, params_as_args, params_as_args2). With no multivariate subclasses left, none of it is reachable: nothing in mpisppy, examples or doc names any of those symbols. So this removes 596 lines that statdist/README.md already said were not here, and the README now says so about base_distribution.py too. The univariate half of the file goes from 42% to 98% covered (what is left is the two abstract method bodies and the interactive plt.show()). Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The epi-spline distributions fit with a nonlinear solver, so with no ipopt on the runner the two epi-spline tests skip and splines.py never runs at all -- Smoothed_boot_epi ships with no verification behind it. This gets ipopt the way Pyomo's own CI does: apt-get the BLAS/LAPACK/gfortran libraries the binary links against, then unpack the IDAES idaes-ext solvers release (30MB) onto the PATH. `ipopt -v` at the end of the step means a bad download fails the job instead of quietly leaving the tests skipped. Verified locally by installing the same tarball: both epi-spline tests pass (serially and under mpiexec -np 2), splines.py goes from 6% to 39% covered and distributions.py from 87% to 97%, since the epi-spline class body now runs. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
UnivariateStudentDistribution named its third argument the variance, and the fit docstring promised "the mean and variance of the distribution as the mean and variance of the data", but the constructor handed sqrt(var) to scipy as the *scale*. A t of scale s has variance s**2 * df/(df-2), so the distribution's variance came out as var * df/(df-2) instead of var -- and under the fit's own df rule of 2v/(v-1) that factor is exactly v, so fitting to data with variance 3.88 produced a t with variance 15.08. The parameter was a scale wearing the name of a variance. The constructor now solves that relation for the scale, so the variance argument is the variance for any df, and df <= 2 is refused: such a t has no finite variance at all, so it cannot be given one. That also settles what the df rule was for. 2v/(v-1) is the df at which a t of *scale one* has variance v, so with the scale derived the fitted scale comes out as one and the rule is self-consistent. What it cannot do is fit a sample variance of one or less -- a t of scale one always has variance above one -- which is what the old "Impossible to define a student distribution" print was about before it went on to build a df=1 Cauchy, a distribution with no variance whatsoever. That case is now a ValueError telling the caller to choose a df themselves or fit something else. Nothing in the repository fits this distribution: the smoothed methods only ever ask for univariate-epispline or univariate-kernel. Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…part 2) Now that the constructor derives scipy's scale from (var, df), any df > 2 reproduces the requested variance, so the old df = 2v/(v-1) is just a free tail-heaviness knob -- and a poor one: it is not scale-free (df changes when the data are rescaled), its coupling runs backwards (larger variance -> df toward 2 -> heavier tails), and it cannot fit a sample variance <= 1. Replace it with method of moments on the excess kurtosis: a t with df > 4 has excess kurtosis 6/(df-4), so df = 4 + 6/excess_kurtosis. This is scale-free and data-driven. A sample with excess kurtosis at or below zero (tails no heavier than the normal) has no finite-variance t that matches it, so df falls back to a large "effectively normal" value instead of raising. The v <= 1 ValueError is gone -- any variance v > 0 fits. Update Test_statdist: drop the 2v/(v-1) df assertion and the low-variance refusal test; add tests that df is recovered from kurtosis (data from a t with df=5 recovers df near 5), that light-tailed data falls back to the normal, and that low-variance data now fits with its moments honored. Nothing in the repository fits this distribution (compute_smoothed_ci only requests epispline/kernel), so the blast radius is confined to these tests. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Second and final PR of the boot-sp move (Stage 1), stacked on #783 (merged).
It adds the dependency-bearing half: the trimmed
statdistlibrary and thesmoothed bootstrap/bagging methods, which completes all eleven
BootMethods.Design:
doc/designs/bootsp_merge_design.md(§2.5, §3, §4, §6).What it adds
bootsp/statdist/— a trimmed port of the univariate distributions(~3,000 of ~5,300 lines). The multivariate/copula/vine/bicop code is left in
the archived repo, which also removes the
from scipy.stats import mvntop-of-file import that scipy ≥ 1.14 breaks. scipy is imported lazily via
pyomo.common.dependencies, so the empirical methods stay scipy-free.smoothed_boot_sp.py— the fiveSmoothed_*methods, fitting aunivariate distribution to the sample and resampling from the fit, plus a
single
compute_smoothed_cidispatch point shared byuser_bootandsimulate_boot. This activates the smoothed branch ofsimulate_bootthatPR-1 shipped behind a "not yet merged" error.
farmer,cvar,multi_knapsack. These could not shipin PR-1 because they instantiate statdist distributions on the empirical
path too, so importing them pulls statdist in.
test_boot_sp_smoothed.py(statdist units, theempirical methods on the statdist-dependent examples, and the smoothed
methods), wired into the
confidence intervalsjob andrun_coverage.bashserially and at
mpiexec -np 2in the same commit;run_all.pyentries; andthe smoothed/statdist sections of
boot_sp.rst.Corrections to the ported algorithms
Reviewed against Chen & Woodruff (2024), Distributions and Bootstrap for
Data-based Stochastic Programming. Three deviations from Algorithms 3 and 5:
smoothed_resample_helperadvanced the record index by one per batch while taking a block of
subsample_sizeconsecutive records, so consecutive batches shared all butone of their draws — a sliding window, not independent resamples.
Algorithm 3 draws a fresh set of
Npoints per batch. Striding by the batchsize (as
smoothed_baggingalready did) widens the interval about fivefoldon the
cvarexample atN=20, nB=10; the previous intervals were far toonarrow.
center_smoothedran whileuse_fittedwas stillFalse, so the center was the raw-data gap.Algorithm 3 takes it from Algorithm 2 run on the fitted distribution, and
that smoothed center is the paper's leading conclusion.
np.varusedddof=0. The algorithm specifies sample variances, andthe empirical estimators already use
ddof=1. At theB_I=3used inrun_all.pythis understateds1by a third.B_I < 2leaves thatvariance undefined and was silently yielding a
nan; it is now an error.Also here: maximization now raises instead of reporting
[0, 0]. Theestimators form every gap as (value at
xhat) − (optimal), which isnon-positive for a maximization, and the drivers floor the interval's lower end
at 0. Per the repo rule that maximization either works or errors,
solve_routinerefuses it, and the message says what real support would take.
The port fixes carried over from the design doc's §4 list are recorded there
(items 13–14 cover the two algorithm corrections above).
The design doc's Status block is refreshed here as well: it still described
PR-1 as an open draft, and now records PR-1 (#783) and its
test_boot_sp.pynp=2follow-up (#820) as merged, with PR-2 being this branch.Validation
ruffclean; the bootstrap suites pass serially (37 passed, 2 skipped) andunder
mpiexec -np 2(22 passed, 2 skipped); docs build with no warnings. Thetwo ipopt skips are the epi-spline fits, which need a nonlinear solver.
No locked values move: the smoothed tests assert structure (ordering,
finiteness, coverage bounds) rather than digits, and PR-1's empirical locked
values are untouched.