Release v0.19.0 - #336
Merged
Merged
Conversation
Weibull.fit on three tied observations raised an IndexError from numdifftools, four steps from the cause. A probability plot has no slope through a single distinct abscissa, so polyfit returned a nan; that nan seeded the maximum likelihood fit, which started at nan and produced a nan hessian; the numerical fallback then asked numdifftools for one, and its list of finite-difference steps came back empty. Neither truncation nor censoring was involved, despite where the symptom was first seen. Gamma and Beta failed the same way but silently. Their moment-based initialisers divide by a variance that is exactly zero for tied data, giving (inf, inf), and a failed optimiser reports its initial guess (#261), so those infinities came back as a fitted model. The probability-plot regression falls back to a unit slope through the centroid when rank deficient. Zero slope is the more literal reading of "no information", but every unpack_rr divides by the slope to recover a scale, so zero only moves the nan one step later. Unit slope leaves each distribution's own unpack_rr in charge of producing correctly typed parameters. Gamma and Beta seed the exponential and uniform cases rather than dividing by zero, and a fit now refuses to return any non-finite parameter whatever produced it. Fixing the nan was not enough on its own: the same fit then returned beta = 512 with success=True and no warning, the optimiser having walked part way along an unbounded direction. So a fit is also rejected when the data cannot pin down the free parameters -- fewer distinct non-right-censored values than free parameters means a flat, or here unbounded, direction in the likelihood. The count is of free parameters, not of the distribution's, so fixing one buys back a degree of freedom: Weibull.fit([10.], fixed={'beta': 2}) is well posed and now returns alpha = 10 where it used to raise. That is why this cannot be a per-distribution constant. One-parameter distributions are untouched: Exponential and Rayleigh fit tied data correctly and must not be caught. Probability plotting is exempt, being a regression rather than a likelihood maximisation, and being how several distributions seed themselves -- that internal call does not carry the caller's fixed, so checking it rejected well posed fits. 330 reference fits across thirteen distributions, five methods, and plain, right-censored and offset data are bit-identical. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Reject degenerate data instead of dying inside numdifftools
Every maximum likelihood fit ran five optimisers -- Nelder-Mead, Powell, BFGS, TNC and Newton-CG -- and kept the best. Over 102 fits across eleven distributions, five data shapes and two sample sizes, all five agreed on the objective to 1e-10, so the last four were confirming what an earlier one had already found. The confirmation was expensive. Nelder-Mead and Powell are derivative free and pay for their robustness in function evaluations -- 50 and 22 against BFGS's 21 -- and each evaluation is O(n). On a million observations those two alone were 42% of the fit. Gradient methods now run first and the loop breaks at the first that converges. Order and early exit had to change together: stopping early without reordering halts at Nelder-Mead, the most expensive rung and the one with the worst objective, while reordering without stopping early saves nothing. The derivative-free pair remains as the fallback and still starts from the cold initial guess, so the multi-start behaviour that motivated the original order survives for the fits that need it -- they are no longer paid for by the fits that do not. Cold-start BFGS wins 83 of the 102 fits, TNC 10, Newton-CG one. The eight that Nelder-Mead or Powell used to win now land on a gradient method at the same objective; they were winning ties on ordering, not finding better optima. 102 of 102 objectives are unchanged, none worse, and the sweep runs 2.20x faster (17.35s against 7.87s). This is not bit-identical. Fitted parameters move in about the seventh significant digit -- median 3e-8, 90th percentile 8e-7 -- because a different optimiser stops at a marginally different point on the same optimum. One offset ExpoWeibull fit moved 13%, on a flat ridge, to a better likelihood (84.5378 -> 84.5276). The GeneralizedOneRenewal docstring and the two tests that pin its simulated MCF are regenerated accordingly. Those values were a snapshot of the library's own output, taken in 00fcfdc, rather than an external reference. Their tolerance is deliberately left tight so that a future change to the optimiser surfaces as a decision rather than passing unnoticed, exactly as it did here. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Stop the MLE ladder at the first optimiser that converges
The beta survival tree and forest tests were 97 of the suite's 180 seconds for 85 of its 2000-odd tests. They now need --run-ml, which continuous integration passes, so the default local run drops to about two minutes with no loss of coverage. The same mechanism carries a new sweep, --run-invariants. Every defect found this cycle slipped past the whole suite, and each lived at an intersection of dimensions the suite tests one at a time: a censored observation that was also truncated, an offset combined with a particular method, an offset combined with a large shift magnitude, and a sample below the 1000 floor of FIT_SIZES. The suite covers each axis; none of it covered the crossings. The full cross is around 600,000 cells, so this does not attempt it. It asserts cheap invariants over a seeded sample instead -- finite parameters, finite neg_ll, a survival function in [0, 1] that never increases, and MLE attaining the lowest negative log-likelihood of the five methods. Four of the five defects would have failed the first two assertions, without anyone having to guess where to look. Accuracy is deliberately not asserted; parameter recovery is test_fit.py's job, and this file is for answers that are not answers at all. A raised exception counts as a pass. Many cells are legitimately refused -- a distribution that cannot be offset, data that cannot identify its free parameters, a method that does not support left censoring -- and only a *returned* model is held to the invariants. Data scale is an axis because nothing tested it before, despite the MLE failure warning advising users to rescale towards 1. It earned its place immediately: MLE at scale 1e3 takes 0.998s against 0.013s at scale 1, independent of n, affecting no other method. That is pre-existing -- verified against 2bacc3b~1, where the same fit took 0.885s -- and is filed for its own investigation rather than fixed here. An earlier draft swept sizes inside each cell, which cost 3,750 fits and never finished inside fifteen minutes. Sizes now have their own test at n = 1, 2, 5 and 40, and the pairwise structural test is pinned to scale 1, since it is about keyword interaction rather than magnitude. 270 cases, three and a half minutes, all passing. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
The previous commit made the beta tree and forest tests skip by default, which would have silently dropped 85 tests from the deployment run. Continuous integration passes --run-ml so its coverage is exactly what it was; only the default local run gets shorter. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
Any fit with a truncation bound on one side only had a nan gradient. ll_interval_or_truncated chose between the CDF and an analytic limit with np.where, which selects the right value but evaluates both branches, so ff(inf) was still taped by autograd and its nan derivative propagated through the selection whichever side won. The objective was always correct; only the gradient was nan. So BFGS and Newton-CG each failed after one evaluation, TNC burned its full 1000 evaluation budget reaching the same conclusion, and Nelder-Mead finished derivative free. A Weibull that fits in 0.014s took 1.36s. A tl of 0 -- a no-op, since F(0) = 0 -- cost exactly as much as a real truncation, which is what gives the cause away, and windows with both bounds finite were always fast because no infinity ever reached the tape. The infinity is substituted out of the argument before the CDF sees it, so one vectorised call covers every row whatever its pattern of bounds, and the surviving np.where only chooses between values that are already finite. xl and xr are data, never traced, so the substitution is invisible to autograd. The stand-in cannot be an arbitrary constant. Zero looks natural and is wrong: a Weibull with beta < 1 has an unbounded density derivative at the origin, so ff(0) would swap one nan gradient for another. Reusing a bound that is genuinely present keeps the stand-in inside the support and at the data's own magnitude. Its value never reaches the result -- np.where discards it -- only its derivative has to be finite. Checked against a truncated Weibull(10, 0.6), which fits cleanly on BFGS. The truncation correction also depends only on the observation window, not on where in it the observation fell, so SurpyvalData now collapses the truncated rows to their distinct windows once at construction. Truncation is nearly always common to the whole sample, so this turned 360 CDF evaluations per likelihood call into one in the test case, and the likelihood runs hundreds of times per fit. Left truncated 1.398s -> 0.022s, right truncated 1.344s -> 0.024s, tl = 0 1.362s -> 0.041s. All 330 reference fits across thirteen distributions and five methods are bit-identical, and BFGS now wins every truncated fit where Nelder-Mead used to. Truncated fits have therefore had no working gradient until now, so cov_matrix and confidence bounds on truncated data were produced from a nan-poisoned hessian path. Whether they were degraded is not established here. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
The sweep is a net for exploring, cast deliberately when the fitting paths are being worked on. Three and a half minutes on every pull request across three Python versions buys little when its assertions hold, so it stays a local target and the conftest no longer claims otherwise. Continuous integration still passes --run-ml. The previous commit speculated that truncated fits, having had no working gradient, might have produced degraded confidence bounds. They did not. mle.py already recomputes a numerical hessian whenever the autograd one is nan or asymmetric (#270), so it caught this on every truncated fit: checked against 906f0cb~1, the standard errors are identical to eight figures. The fallback was part of what made those fits slow, which the gradient fix has now removed. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
Every parameter bounded on (0, inf) is mapped to the unbounded space the optimiser searches by adj_relu, which chose between x + 1 and exp(x) with np.where. Autograd evaluates both branches, so exp(x) was taped even where x + 1 was selected, and above x = 709.78 it overflows to inf -- poisoning the derivative of the branch that *was* chosen. The transform's jacobian came back nan, and cov_matrix is that jacobian either side of the inverse hessian. A Weibull fit to lifetimes in hours had standard errors; the same lifetimes in seconds returned nan for all of them, with no warning. The threshold is a property of the fitted parameter, not the sample size or the conditioning, which is why it looked arbitrary. A fit died as soon as any (0, inf) parameter exceeded about 710, and that predicts all twelve distributions exactly: Weibull, LogLogistic and ExpoWeibull at scale 100, Rayleigh at 141x, the Gumbels and Logistic at 350x, and the Gamma at the *small* end, its rate growing as the data shrinks. Normal, LogNormal and Exponential were immune because they have closed-form estimators and never touch the transform; the Uniform reports no covariance at any scale by design, its MLE being an order statistic. The hessian was never the problem. The numerical fallback (#270) returned a finite, well conditioned matrix at every scale -- condition number 2.41 for the Gumbel and 1.0 for the Rayleigh at 1e6 -- which is why this presented as silent nan rather than a warning or a failure. It also explains the speed: the same nan reached the objective's gradient, so BFGS, TNC and Newton-CG each gave up and Nelder-Mead finished the fit derivative free. Twelve distributions at seven scales go from 19.59s to 1.18s. Restoring the gradient exposed a second scale problem underneath. scipy stops BFGS when max|grad| < gtol, and its default of 1e-5 is an absolute threshold on a quantity that is not scale free: a log-likelihood's gradient shrinks like 1/theta, so on data in the tens of thousands the test is met well short of the optimum. This was invisible while those fits ended on Nelder-Mead, which is derivative free and kept going. Three reference fits landed 1e-2 away, at a likelihood 2e-3 below the answer they had been recorded from. gtol is now scaled by the gradient where the search starts, which asks that it fall by a fixed number of orders of magnitude -- the same requirement at every scale. A fixed tighter constant is not equivalent: 1e-8 fixes the large fits but is unreachable for small samples, dropping an n=8 Weibull into TNC. Results that already worked are unchanged: 65 of 96 reference fits are bit-identical and the other 31 are the restorations, agreeing on the objective to fifteen significant figures. The restored standard errors satisfy se(theta * c) = c * se(theta) to 1e-5 or better everywhere tested, most to 1e-8, and to machine precision for the Rayleigh and Normal. Fits now converge slightly further wherever BFGS wins, always towards a better likelihood, so a few pinned numbers moved in their last digits. The two Monte Carlo tests in test_counting.py are loosened to rtol=1e-3: they drive a 5000-run simulation from an optimiser's output, where a seventh significant figure moves the MCF in the fourth. Also moves conftest.py to the repository root. pytest_addoption is only honoured in initial conftest files, and the CI invocation selects with --ignore and names no path, so the flags were never registered and the run would have died on "unrecognized arguments: --run-ml". #323 is rescoped rather than closed. It had proposed internal rescaling to fix both the bounds and the speed; neither needs it. Re-measured with the gradient working, rescaling is between 4.2x faster and 6x slower depending on the distribution. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Fix two autograd np.where traps: truncated fits and the parameter transform
The changelog said the restored standard errors satisfy se(theta * c) = c * se(theta) "to 1e-5 or better across every distribution and scale tested". That was measured against a fixed gtol of 1e-8, which was then replaced by a relative one before the sweep was re-run. It does not hold for what shipped: the Weibull at data scale 1e6 sits at around 1e-3. Re-measured against the merged code, the true figure is better than claimed almost everywhere -- 1e-9 or better up to scale 1e4, several distributions at 1e-11, machine precision for the Rayleigh and Normal -- with a single outlier that the old wording quietly covered up. The outlier is the interesting part, so it is now stated rather than averaged away. Scaling gtol by the gradient at the initial guess does not make the criterion scale free the way it appears to: the initialiser scales with the data too, so the gradient at the starting point is itself roughly scale invariant. A Weibull at scale 1, 1e4 and 1e6 all come out with gtol = 1.86e-6. It bought a better constant, not scale independence. Nor is there a constant that serves every case. Measured across the three large-magnitude reference fits, an n=8 sample and a Weibull at 1e6: 1e-8 relative fixes the last but drops the small fit into TNC 5e-5 away, and 1e-9 pushes both small fits to TNC. 1e-7 is simply where everything the suite covers passes. Recorded in #323, which is now scoped to preconditioning the optimiser's search -- the cheap form, touching neither the reported parameters nor the covariance, and so needing none of the twelve per-distribution back-transforms that made the original proposal risky. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Correct an overstated claim about scale-equivariant standard errors
scipy stops BFGS when max|grad| < gtol, an absolute threshold on a
quantity that is not dimensionless. A log-likelihood's gradient moves
with two things at once: it shrinks like 1/theta with the units the data
is recorded in, and it grows like n, being a sum over observations. The
default of 1e-5 therefore means something different for every dataset,
and both directions bite.
On data measured in tens of thousands the test is met well short of the
optimum. Three reference fits landed 1e-2 away in relative terms, at a
likelihood 2e-3 below the answer they had been recorded from -- hidden
only because those fits used to end on Nelder-Mead, which is derivative
free and kept going. At n = 1e5 the gradient is five orders of magnitude
larger, the same threshold is unreachable, and BFGS gives up on censored
samples it should handle easily. The second direction was found by
benchmarking against lifelines with censoring, where it was the single
configuration in 64 where surpyval was slower.
Tuning the constant cannot fix a criterion whose meaning moves. Four
alternatives were measured against the same reference set before
settling on this. Scaling gtol by the gradient at the initial guess
looks scale free and is not -- the initialiser scales with the data too,
so that gradient is itself roughly scale invariant, and a Weibull at
scale 1, 1e4 and 1e6 all came out with gtol = 1.86e-6. BFGS's xrtol and
L-BFGS-B's ftol are scale free in form but stop on how the optimiser is
behaving rather than on the quantity that is zero at the answer, so they
quit early along flat directions -- precisely where the standard error
is largest and most needs to be right. The ExpoWeibull, three parameters
and a flat surface, was 17% out under both.
So rescale the problem in both of its dimensions and let a single
constant mean the same thing everywhere:
s = max(|u0|, 1) f0 = max(|f(u0)|, 1)
v = u / s g(v) = f(s v) / f0
The starting point is order 1 in every component and so is the
objective, whatever the units and whatever the sample size. Since
dg/dv = s * df/du / f0, with s growing exactly as df/du shrinks and f0
growing like n, so is the gradient.
scale-equivariance of se relative gtol rescaled problem
worst deviation 1.6e-3 2.0e-5
n=1e5, 30% cens, scale 1 0.163s 0.153s
n=1e5, 60% cens, scale 1e6 0.250s 0.119s
Rescaling the parameters alone reaches 2.8e-8 on the first row but is
the version that leaves BFGS failing at large n, where those two fits
take 0.370s and 0.590s. Normalising the objective as well trades a
little of that accuracy for a criterion that holds across sample sizes.
Both mappings are fixed before the search begins and neither can move
the optimum: a diagonal linear change of variable relocates a minimum no
more than dividing the objective by a positive constant does. res.x and
res.fun are both mapped back inside the helper -- res.fun because the
ladder compares it across rungs -- so no scaled quantity exists anywhere
else in the package, not even transiently. The covariance step builds
its own hessian at the returned point and sees exactly what it saw
before; scipy's res.hess_inv would be in scaled units and is read
nowhere.
One known wrong answer comes with this. WeibullPH.fit_tvc returns a
degenerate fit on the fixture in test_episode_split_is_an_identity,
marked xfail(strict=True) against #326 so it reverses on its own when
fixed. Nothing here touched the regression path, which has its own
optimiser; the seed it receives moved by six parts in ten million and
that was enough. Its first stage still finds the right answer (neg_ll
927.83) before TNC walks into the region where a left-truncated
likelihood is unbounded, returns -21311.40, reports success, and is kept
unconditionally by optimise_ph. Neither a success check nor a comparison
of objectives would catch it, the divergent answer having the lower one.
That fit has been beside this basin all along with only its starting
point keeping it out.
Full suite plus the invariant sweep: 2368 passed, 6 xfailed.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…flows
A parametric regression fit to left-truncated data could report a
log-likelihood tens of thousands higher than its parameters earn, and be
optimised towards it. A WeibullPH fit reported neg_ll -21311.40 at
parameters whose true value is +17118.30, against 927.83 at the correct
answer.
truncation_correction computed the mass in each observation's truncation
window as a difference of CDFs, floored at the smallest positive float:
np.log(np.maximum(right - left, _TINY))
Under left truncation that difference is 1 - F(tl), the survival
probability at the truncation bound, which underflows to exactly zero as
the fitted scale shrinks. The floor then capped the correction at
log(tiny) = -708 instead of letting it grow without bound. The caller
subtracts the correction, so every truncated row appeared to contribute
+708 to the log-likelihood: a region the data rules out entirely became
the best fit on offer.
Evaluate one-sided windows in log space through log_sf and log_ff, which
stay finite where the difference cannot, so there is nothing to floor.
Only a genuinely two-sided window still takes a difference, and there
both bounds are finite and the mass is not driven to zero by the scale
alone. The np.where branches evaluate at substituted-finite arguments so
an infinity in an unselected branch cannot poison the gradient of the
selected one -- the same care the univariate likelihood needs. _ff_safe
existed only for the old formulation and is removed.
This lives in _likelihood.py, shared by proportional hazards,
proportional odds, accelerated failure time and accelerated life, so any
left- or right-truncated parametric regression fit was exposed. It
needed only the optimiser to wander far enough for the underflow to
bite, and nothing warned when it did.
The first diagnosis of this was wrong and is worth recording. It looked
like a genuinely unbounded truncated likelihood being followed
legitimately by an optimiser, which would have made it a statistical
pathology needing bounds or rejection. Recomputing the likelihood by
hand from the model definition disproved that: at the correct parameters
surpyval agrees to the digit, so the objective is right everywhere the
floor does not engage, and the degenerate point is simply much worse
rather than much better. The bug was arithmetic, not statistics.
The xfail added against #326 in the previous commit is removed, and the
regression test pins the defect directly rather than relying on the
identity test to catch it: the reported objective must equal an
independently computed one, and shrinking the scale ten, a hundred and a
thousandfold below the truncation bounds must always score worse.
Full suite plus the invariant sweep: 2370 passed.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
An entry time below every observation excludes nobody, so it should
leave a fit untouched. With any left-censored row present it raised
instead:
ValueError: An observation's censoring interval does not intersect
its own truncation window ...
A support index j stands for the half-open interval
(bounds[j], bounds[j+1]], so an event placed there is already strictly
after bounds[j]. The first index a row entering at tl may use is
therefore the last bound equal to tl -- that interval is (tl, next].
The window construction took one index further on and discarded it.
This mattered most for left censoring because such an event lies in
(-inf, xr], which under an entry at tl is the single interval (tl, xr],
frequently the only one the row has. Dropping it left the row with an
empty support.
Neither endpoint of the search alone is correct, which is what made the
off-by-one awkward to see. side="left" keeps the zero-width (tl, tl]
interval that a duplicated exact event time creates, readmitting an
event at exactly the entry time and breaking the strict (entry, exit]
convention (#260) -- three tests catch that. side="right" discards
(tl, next] as well, one interval too many. side="right" - 1 lands
between them, and exactly, because every finite truncation time is
itself in bounds.
Half of #308 remains. Convergence on data mixing all four censoring
types with left truncation is unaffected, and the cause is not the one
the issue supposed: expected event counts come back inflated, about 3.1
events at every observation time for thirty observations, so the
estimator ladder exhausts its risk set within a handful of steps and the
survival curve collapses. That is the ghost/normalisation step, not the
index arithmetic.
The new part is the trigger. A common entry time round-trips exactly;
entry times that differ are what activate it, since only then does a
later-entering row have unseen deaths imputed below its own window. That
narrows the reproducer from thirty points to six, added as a strict
xfail carrying the diagnosis. The sweep case keeps its xfail with the
reason narrowed from "left censoring + left truncation" to convergence
alone.
Full suite plus the invariant sweep: 2373 passed, 5 xfailed -- the same
five as before, now pinned to a narrower defect.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
Left-censored observations combined with two or more distinct entry times admit a flat direction in the likelihood, and where the data leans on it the estimate is worthless while looking ordinary. A six-point example returns sf(2.5) = 0.002 where the sensible answer is 0.78. An interval that one observation could have failed in, but that precedes another observation's entry, is worth mass to the first and costs the second nothing: the second's contribution is conditional on its own entry, so mass it never had the chance to see divides out of both its numerator and its denominator exactly. Its numerator falls to 1.5e-5 and its denominator to 5e-5, and the ratio is unchanged. The estimator drives 99.995% of the mass into a single such interval, reaching a log-likelihood of -6.14 against -9.36 for the sensible answer. So the EM is not misbehaving. It maximises correctly and the likelihood has no interior maximum to find, which is why raising max_iter never helps. Both ingredients are needed: left censoring, the only kind whose support reaches back into the entry region, and two distinct entry times, so the interval exists at all. Six distinct entry times without left censoring fit flawlessly; one common entry time with left censoring round-trips exactly. The fit is still returned. Meeting the condition does not mean the data is spoilt -- across 240 simulated samples that all met it, the proportion actually degenerating ran from 8% to 72% with the share left censored, so rejecting on structure would refuse far more good data than bad. What separates them is how much mass lands on the flat direction: healthy fits reached at most 0.836 of it, spoilt ones a median of 0.994. Warning above 0.9 caught them with no false alarms across those samples, where 0.7 would have cost 9% and 0.5 40%. The share is exposed as model.exploitable_mass so a borderline fit can be judged rather than guessed at. The structural condition is deliberately not a trigger by itself. Ordinary staggered-entry data meets it routinely and estimates perfectly well -- an event earlier than a later unit's entry also sits outside that unit's window, but it is pinned by the data rather than free to move. An earlier draft fired on the flag plus non-convergence and would have told the #203 case that more iterations would not help, which is false: that case is structurally exploitable and converges once given them. Two existing tests caught it. Full suite plus the invariant sweep: 2378 passed. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
The three remaining xfails are the non-identifiable configurations: left censoring plus two or more distinct entry times, where the likelihood has a flat direction whose supremum sits on the boundary of the parameter space. That is not a bug in the fitter, so #308 closes with the warning that now reports it; the question of what such a fit should *return* moves to #327, where the Vardi/Wang connectivity criterion is the candidate replacement for the mass heuristic. Also adds coverage for the warning itself: that a genuinely non-identifiable fit is reported, that an identifiable one is not, and that distinct entry times alone are harmless. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
The PH objective closes over regression_neg_ll, which is written in autograd.numpy and is therefore differentiable — but the ladder never passed a jac, so scipy fell back to a two-point finite difference and paid p + 1 extra objective evaluations per gradient. Nor was the search preconditioned, so PH kept the scale sensitivity the univariate MLE ladder was cured of: at data scale 1e6 a WeibullPH fit settled 1.5 nats of log-likelihood short of the optimum and 1e-2 away in the covariate coefficients, which is a different fitted model rather than a tolerance artefact. And TNC's answer was returned whether or not it had converged, so a rung that can only ever be an improvement was free to be a regression — the AFT/PO ladder defined immediately below already guarded against that. The ladder is now preconditioned BFGS on the analytic gradient, then TNC, then Nelder-Mead, stopping at the first rung that converges and never returning a point worse than it started from. Against lifelines, scoring both answers on an independently written Weibull PH log-likelihood, the scale-1e6 shortfall goes from 1.468 nats to 1.1e-08 and a 50000 x 10 fit drops from 1.53s to 0.40s (lifelines 2.38s). preconditioned_bfgs moves from fitters/mle.py up to fitters/__init__.py beside bounds_convert and fallback_minimize so both ladders share one copy. The univariate ladder is unchanged. optimise_nm_tnc, which serves AFT and PO, has the same two omissions. Its Nelder-Mead first rung means the failure mode may differ, so it is left to be measured rather than changed blind. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
Four costs, none of them the maths — the coefficients still match lifelines to between 1e-06 and 1e-12, digit for digit with before. _GroupBy.sum used np.add.at for its multi-dimensional case, an unbuffered scatter with no fast path, called about ten times per jac_hess on (n, p, p) arrays: 6.3s of a 16.9s fit at n=50000, p=10. It is now a sorted np.add.reduceat, and skips even that when the keys arrive already grouped or all distinct — one-element groups in order are the input array, which is what continuous event times give you. The hessian was a Python double loop with an np.outer inside it, 4.8s of the same fit plus 370000 calls to np.outer. The sum over tied deaths factors out of the p x p part: only c = j/d depends on j, so five scalar sums per event time carry the ragged axis and the covariate blocks are formed once. O(times x ties x p^2) becomes O(times x ties + times x p^2), and the no-ties case stops needing a branch — an untied time is a single j = 0 term with c = 0. Z's outer product was rebuilt inside jac_hess every iteration though Z is fixed for the fit; it is hoisted, and the two weighted-outer einsums are plain broadcasts. Rows are put in event-time order once when the closures are built so _GroupBy never permutes an (n, p, p) array. Nothing downstream depends on row order and the model still stores the caller's unsorted arrays. Efron, 40% censored: n=50000 p=2 from 5.71s to 0.66s (lifelines 3.13s), n=10000 p=5 from 1.38s to 0.34s (0.66s), n=2000 p=10 from 0.38s to 0.11s (0.16s). Heavy ties gain most: n=20000 p=5 from 1.91s to 0.33s. Delayed entry at 43384 rows from 6.86s to 2.62s, and 20000 start-stop rows from 1.07s to 0.44s. Two cases stay slower than lifelines — n=50000 with p=10 (8.3s against 4.0s) and heavy ties (0.33s against 0.08s). Both are now dominated by materialising (n, p, p) arrays, 40MB apiece and several per iteration. Getting past that means accumulating the information matrix per event time rather than building per-observation outer products, which is a change of algorithm and belongs in its own change. The algebraic collapse is pinned by oracle tests that run the loop it replaced, over no ties, heavy ties, empty times, and the fractional count weights where range(int(d)) and c = j/d disagree. Also closes a literal block in the preconditioning changelog entry that was rendering |u0| as an undefined substitution. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Truncated-likelihood and Turnbull correctness fixes, and a faster PH stack
Both were written before the issues existed, so they described work that was deferred without saying where it had gone. The Cox (n, p, p) materialisation is #332; measuring optimise_nm_tnc before changing it is #331. The AFT/PO sentence also now says *why* it is being measured rather than fixed by analogy: Nelder-Mead is derivative free, so it cannot fail the way BFGS did on large-magnitude data, and assuming it needs the same change would be a guess. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Point the changelog's two follow-ups at their issues
The batch is not patch-shaped. Fitted parameters move — the truncated regression fix was worth 38430 log-likelihood on the case that exposed it, and preconditioning changes any fit on large-magnitude data, so pinned outputs will differ. Turnbull now fits left-censored data under left truncation where it used to raise, emits an identifiability warning on data that previously fitted silently, and returns a new exploitable_mass field. Degenerate data raises a different exception than it did. None of that is a new capability anyone asked for, which is the case for calling it a patch. But it all changes underfoot: same code, same data, different answer. The repo's own precedent agrees — 0.15.1 was a single hang fix and 0.15.2 four small hardening items, while 0.18.0 carried a comparable batch of fixes and perf work as a minor. SCHEMA_VERSION in serialisation.py is deliberately untouched: it tracks the serialised model format, which has not changed, and is not coupled to the package version. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Number the pending release 0.19.0, not 0.18.1
Date the changelog heading. The version itself was set in the previous commit, along with the reasoning for a minor rather than a patch. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
…uatqe Cut v0.19.0
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
2 participants
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.
Everything on
developsince v0.18.0.Correctness
Confidence bounds no longer turn to nan on data in large units (#4ce0d1d). Every
(0, inf)parameter goes throughadj_relu, which chose betweenx + 1andexp(x)withnp.where. Autograd evaluates both branches, so abovex = 709.78the unselectedexpoverflowed to inf and poisoned the derivative of the branch that was chosen. A fit died as soon as any such parameter exceeded ~710 — which is why the threshold looked arbitrary. Also ~17x faster, because the same nan reached the objective's gradient and every optimiser in the ladder gave up.The truncation correction could manufacture likelihood out of an underflow (#326). The CDF difference was floored at
_TINY; under left truncation that difference isS(tl), which underflows to zero on a wide window, capping the correction atlog(tiny) = -708— and the caller subtracts it, so every truncated row gained 708 log-likelihood. Worth 38 430 on the case that exposed it.Turnbull kept the interval starting at an observation's own entry time (#308).
searchsorted(..., side="left")is one atom too high whentllands on a support boundary, so a left-censored row entering at an event time could not have failed in the interval beginning there.A Turnbull fit resting on a non-identifiable direction now says so (#308). With left censoring and two or more distinct entry times the likelihood has a flat direction whose supremum sits on the boundary — the degenerate answer scores −6.14 against −9.36 for the sensible one, so this is the likelihood's behaviour, not a fitter bug. The fit warns and reports
exploitable_mass. The trigger is mass evidence only: over 240 fits, only 8–72% of structurally eligible configurations actually degenerate, so rejecting on structure would refuse ~88% of good data.Parametric PH degraded on large-magnitude data (#328).
optimise_phpassed nojacdespite the objective being autograd-traceable, was not preconditioned, and returned TNC's answer whether or not it converged. At data scale 1e6 aWeibullPHfit settled 1.5 nats short and 1e-2 away in the coefficients.Degenerate data is rejected with an explanation instead of dying inside numdifftools.
Performance
1/thetaand grows liken, so no absolutegtolserves every fitAgainst
lifelines, scoring both packages' answers on independently written log-likelihoods: surpyval now wins on delayed entry (2.62s vs 170s at 43 384 rows), time-varying covariates (0.44s vs 7.8s at 20 000 rows), and every parametric PH cell measured.Known limitations, tracked
lifelines(n=50 000 with ten covariates, and heavy ties). Both are dominated by materialising(n, p, p)arrays; getting past it means accumulating the information matrix per event time, which is a change of algorithm.Five xfails in the suite are the non-identifiable Turnbull configurations, marked strict so that a future change making one converge sensibly will be reported.
Why minor rather than patch
Every item is a fix or a speed-up, but fitted parameters move, Turnbull accepts input that used to raise and emits a warning on data that previously fitted silently, there is a new output field, and degenerate data raises a different exception. The changes are invisible at the call site: same code, same data, different answer.
SCHEMA_VERSIONstays at 1 — the serialised model format has not changed.Verification
Full suite 2393 passed, 1 skipped, 5 xfailed. CI green on every commit: lint (flake8, mypy, black) and pytest on Python 3.11, 3.12 and 3.13.
🤖 Generated with Claude Code
https://claude.ai/code/session_01TRhKL2fJBiAAfNo9rihwts
Generated by Claude Code