fix(engine): correct BoxCox lambda selection (diffuse Jacobian + vari…#13
Merged
Conversation
…ance ridge) The joint/profile BoxCox lambda MLE was biased, overshooting toward a bound (iid N(1000,50) -> lambda=2 not 1; Nile -> lambda=2 + spurious global trend). Two coupled defects in the concentrated likelihood (llik, SSpace.h): 1. Diffuse observations received their full-scale BoxCox change-of-variable Jacobian (lambda-1)*sum log(y) but NOT the concentrated-variance scale term -0.5*nDiff*log(s2) (dropped earlier to dodge a +Inf logLik). The KF runs in concentrated units, so log|Finf,t| = log(s2) + log(Finf*,t); dropping log(s2) left the Jacobian dangling -- a monotone-in-lambda term growing with the number of diffuse states, biasing lambda. Restored it: this is the exact-diffuse concentrated BCnorm marginal likelihood. 2. Restoring log(s2) reintroduced the s2->0 singularity on degenerate fits (short over-parameterised series interpolating the data -> +Inf logLik). Replaced the former hard floor with a weakly-informative inverse-gamma prior on the concentrated variance: s2 = (SSR + nu0*s0sq)/(n + nu0), nu0=1, s0sq = 1e-4*Var(g(y)). Bounds s2 smoothly; healthy fits (SSR >> nu0*s0sq) untouched; AICc's k-penalty + the existing n-k-1>0 admissibility do the rest. Verified: lambda recovers truth on clean oracles (iid->1, mult->0, additive->1) and real data (Nile-> local level lambda=1, AirPassengers-> log lambda=0, UKgas-> ~log); short-series singularity resolved (logLik negative); selections stable across prior strength 1e-6..1e-4. Shared engine, so identical in R/Python. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh
…ence bands
The engine computed component state variances internally but never returned
them ("dead compV"). Surface them through the R interface so component
uncertainty bands (and the prediction-variance widening over missing blocks)
can be plotted.
- components() gains a smoothedVar flag. Default (false) keeps the fast state
smoother -> filtered variances (cheap; enough for the widening-over-a-gap
picture). When true, the full backward Nt recursion runs so inputs.P holds
the two-sided smoothed variances P_{t|T} (O(m^3), gated behind the flag).
- compV threaded MuseOutputs -> musecpp2R -> pts object as $compV (BC scale,
columns Level/Slope/Seasonal/Irregular), reshaped like $comp.
- New pts() argument componentVariance = FALSE; TRUE requests smoothed bands.
- One appended .UCompC arg (compVarSmootheds); RcppExports regenerated.
Verified: filtered/smoothed compV finite & non-negative; smoothed < filtered
(more information); filtered Level variance widens over an artificial gap and
re-narrows. Lambda fix from the previous commit unaffected.
Note: R front-end only -- musecpp2py.cpp / pts.py not yet mirrored (engine
change is backward-compatible: compVarSmoothed defaults false, so Python builds
and runs unchanged, just without the compV output).
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh
polyStationary/invPolyStationary now return/consume the STANDARD Box-Jenkins AR coefficients (pacfToAr convention, no internal negation), so the engine's internal ARMA parameter equals the reported coefficient -- AR is (1 - phi B), MA is (1 + theta B), both matching stats::arima. Removes the scattered compensating flips that this convention mismatch required: - armaMatrices/armaMatricesTrue: T = +col (AR), R = +col (MA). - musecore.h reporting + parameterValues: no AR sign flip. - setEstimatedParams: no AR *= -1 on the true-param path. - AR seed: sample PACF maps straight into estimation space (no negate). This fixes the AR-on-irregular basin bug: the old seed fed a standard (1-phi B) PACF into the internal (1+phi B) invPolyStationary, stranding a correct +phi seed in the -phi basin. For pure-structural models the optimiser escaped it, but under the (correct) exact-diffuse log(s2) term that basin is a local minimum it cannot leave -- the RW level stays active and the AR flips sign. With the seed in the right basin the level collapses on stationary data and AR recovers +phi. MA seed: clamp the arToPacf input to the invertible region (|PACF| < 1) as the AR seed already does; otherwise a 2-block SARMA MA seed divides by (1 - phi^2) at the unit boundary and crashes invPolyStationary. Verified: AR(1)=0.7 -> +0.67, MA(1)=0.5 -> +0.4x, ARMA(2,2)/SARMA signs correct (match arima); all 96 R testthat pass; lambda oracles correct (iid->1, Nile-> local level, Air->log). CLAUDE.md updated with the full rationale. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh
…am refs
- musecpp2py.cpp + core/pts.py: expose component state variances (compV) and the
component_variance flag, mirroring the R $compV feature (filtered by default;
smoothed P_{t|T} when requested). Same shared C++ engine, so identical output.
- core/pts.py: mirror R's joint-lambda + outlier workaround -- when outliers are
requested with an auto-lambda ("Z") power, pin lambda from a quick no-outlier
fit first, sidestepping the engine's dummy-injection dimensionality bug.
- tests/dump_{pts,outlier}_reference.R: dump nparam(m) (scalar total) instead of
the 2x5 $nParam matrix, so the parity harness compares like-for-like.
Verified: 45 Python functionality tests pass; all six R<->Python parity suites
match to <= 1e-6 (engine, pts, forecast, summary, outliers, diagnostics).
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh
The componentVariance argument added with the compV feature had no roxygen @param entry, so R CMD check flagged "Undocumented arguments in documentation object 'pts'" -> WARNING (CI failure). Add the @param and regenerate man/pts.Rd. Verified locally: R CMD check documentation sections all OK; vignette builds; Python ruff + mypy + pytest (45) all pass. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com> Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh
config-i1
added a commit
that referenced
this pull request
Jul 23, 2026
…e-jacobian fix(engine): correct BoxCox lambda selection (diffuse Jacobian + vari…
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.
…ance ridge)
The joint/profile BoxCox lambda MLE was biased, overshooting toward a bound (iid N(1000,50) -> lambda=2 not 1; Nile -> lambda=2 + spurious global trend). Two coupled defects in the concentrated likelihood (llik, SSpace.h):
Diffuse observations received their full-scale BoxCox change-of-variable Jacobian (lambda-1)sum log(y) but NOT the concentrated-variance scale term -0.5nDifflog(s2) (dropped earlier to dodge a +Inf logLik). The KF runs in concentrated units, so log|Finf,t| = log(s2) + log(Finf,t); dropping log(s2) left the Jacobian dangling -- a monotone-in-lambda term growing with the number of diffuse states, biasing lambda. Restored it: this is the exact-diffuse concentrated BCnorm marginal likelihood.
Restoring log(s2) reintroduced the s2->0 singularity on degenerate fits (short over-parameterised series interpolating the data -> +Inf logLik). Replaced the former hard floor with a weakly-informative inverse-gamma prior on the concentrated variance: s2 = (SSR + nu0s0sq)/(n + nu0), nu0=1, s0sq = 1e-4Var(g(y)). Bounds s2 smoothly; healthy fits (SSR >> nu0*s0sq) untouched; AICc's k-penalty + the existing n-k-1>0 admissibility do the rest.
Verified: lambda recovers truth on clean oracles (iid->1, mult->0, additive->1) and real data (Nile-> local level lambda=1, AirPassengers-> log lambda=0, UKgas-> ~log); short-series singularity resolved (logLik negative); selections stable across prior strength 1e-6..1e-4. Shared engine, so identical in R/Python.
Claude-Session: https://claude.ai/code/session_017YTrDQdw1PnrFfpD5uR9Qh