Skip to content

R Python differences

Jump to bottom
Ivan Svetunkov edited this page May 29, 2026 · 10 revisions

R vs Python differences

This page documents where the R smooth package and the Python port produce bit-identical results and where they don't. It is the running record of the cross-language status as the two implementations converge.

What is identical (machine precision)

The following are bit-identical or agree to within a handful of ULPs between R and Python on the same input series:

Quantity Across all scenarios
coef(m) / m.coef — fitted coefficient vector ✓
m$lossValue / m.loss_value — cost-function value ✓
m$logLik / m.loglik ✓
m$fitted / m.fitted — fitted values ✓
m$residuals / m.residuals ✓
forecast(m, h=…) / m.predict(h=…) — point forecasts ✓
sigma(m) / m.sigma — empirical residual std (R formula) ✓
AIC / AICc / BIC / BICc ✓
vcov(m), confint(m), summary(m) — OM, OMG, and ADAM with initial="backcasting" / "complete" ✓
vcov(m) / confint(m) — ADAM with initial="optimal" / "two-stage" ≤2% — see below
coefbootstrap() — structural parity (same coef cardinality, PSD vcov) ✓
multicov(type="analytical") ≤5%
multicov(type="empirical") — calls the same C++ ferrors backend on both sides ✓

What had to be aligned

The shared C++ core and the matching call sites:

  1. Shared finite-difference Hessian — src/headers/hessianCore.h is the single source of truth for the FD Hessian used by both R (vcov.adam / vcov.om / vcov.omg) and Python (numerical_hessian). Per-parameter relative step hᵢ = ε^(1/4) · max(|xᵢ|, 1) so large-magnitude parameters (initial level / trend / seasonal in B when initial="optimal" or "two-stage") get a meaningful perturbation.

  2. bounds="none" during FI computation. Both implementations disable the cost-function bounds penalty during the FD Hessian so perturbations at boundary parameters don't trip the 1e+300 infeasibility return.

  3. abs(diag(...)) of vcov. Mirrors R/adam.R:5226 — both sides agree on the sign convention for non-PSD FIs.

  4. sigma() formula. Python ADAM.sigma computes sqrt(SS / df) exactly per R's sigma.adam (R/adam.R:4625-4658), with a distribution-specific SS and df = nobs - nparam (Python's nparam already excludes the scale parameter that R counts and then subtracts). The internal optimisation scale R calls m$scale is exposed as m.scale (separate from m.sigma).

  5. OM.sigma / sigma.om. OM has no scale parameter on the probability axis; both languages report sqrt(mean(residuals²)) on the link-transformed scale — formula from oes_old / oesg_old.

  6. dgnorm shape lower bound. Python's optimiser lower bound on the generalised-normal shape is 1e-10, matching R's adam_checkOptimizer. A tighter bound (0.25) was previously used and changed the NLopt Nelder-Mead simplex behaviour, producing different local minima on the same cost surface.

  7. two-stage ≡ optimal for FI. Two-stage produces the same B shape as optimal; both R (R/adam.R near the initialTypeFI block) and Python treat them identically when computing the Hessian.

The one remaining numerical gap

The ULP-level gap on ADAM with initial="optimal" (or "two-stage") on multi-parameter ETS models comes from the linear-least-squares solver used inside msdecompose's default global smoother — and nowhere else.

  • R uses .lm.fit → LINPACK Householder QR (dqrls.f).
  • Python uses numpy.linalg.lstsq → LAPACK dgelsd (SVD-based).

Both solve the same least-squares problem; they round at different bits during elimination/back-substitution. Result: msdecompose's level (and to a much smaller extent trend) differs by ~1 ULP on the data scale (level ≈ 100 for AirPassengers → ~1.4e-14 absolute difference).

That ULP enters the x0 passed to NLopt only when level / trend are in B — i.e. for initial="optimal" and "two-stage". The backcasting and complete paths derive these inside the C++ kernel from the data each call, so the LSQ-solver ULP never reaches the optimiser. OM and OMG likewise don't carry initial states in B. The optimiser is deterministic given the same x0, so this single ULP is the entire source of downstream divergence in vcov / confint / summary for optimal / two-stage.

The downstream amplification chain — LSQ ULP → x0 ULP → final-B ULP → FI numerator → inverse-FI — used to compound this into ~5e-3 diffs on the worst vcov entry. With the per-parameter relative Hessian step (point 1 above) the FI amplification is now small enough that the gap collapses to roughly 2% relative on vcov / SE / confint for the optimal-initials path.

This is a property of the BLAS/LAPACK implementation, not an algorithmic bug. Two different LSQ solvers will always round differently. Closing the residual gap would require either (a) using a closed-form formula for the X = [1, t] case inside msdecompose (intercept + time — trivially bit-identical across implementations) or (b) routing msdecompose's smoother through a shared C++ QR via the same Rcpp / pybind11 bridge used for the Hessian.

Simulation functions and RNG-induced differences

This section covers the sim_* family (sim_es, sim_ssarima, sim_sma, sim_oes, sim_ces, sim_gum) and the simulate() methods (ADAM.simulate, OM.simulate, OMG.simulate).

Default seeds do NOT cross languages

R and Python use different random number generators:

  • R: the Mersenne Twister (set.seed(n)) for rnorm, rt, rgnorm, etc.
  • Python: NumPy's PCG64 default (np.random.default_rng(n)).

Even with identical seeds, the two languages produce completely different draws. The same sim_es(model="ANN", seed=42) call on either side gives an entirely different series. The same applies to ADAM.simulate(seed=42), OM.simulate(seed=42), and so on.

What is preserved across the R/Python boundary is the state-space machinery that consumes those draws — the matrix prep, the C++ adamCore::simulate kernel call, the post-processing (link functions for OM/OMG, scaling for ETS, etc.). These are bit- equivalent on both sides. The only divergence is the random draws the kernel ingests.

The "plug-in numbers" trick

To prove the state-space machinery itself matches across languages, all of the sim_* functions and the simulate() methods accept a callable in place of the distribution name:

# Python
errors = np.linspace(-1.5, 1.5, 60)
def feed(n):
    return errors[:n]
sim = sim_es(model="ANN", obs=60, randomizer=feed, ...)
# R — same vector
errors <- seq(-1.5, 1.5, length.out=60)
ourFunction <- function(n, ...) errors[1:n]
sim <- sim.es(model="ANN", obs=60, randomizer="ourFunction", dummy=1)

Both calls feed the same errors into the same C++ kernel, so the resulting series, states, and persistence vectors agree to floating-point. The r_parity-marked tests in the Python repo (test_sim_es_r_parity.py, test_sim_ssarima_r_parity.py, test_sim_sma_r_parity.py, test_sim_oes_r_parity.py, test_sim_ces_r_parity.py, test_sim_gum_r_parity.py) all use this trick and pass at atol=1e-9 or tighter:

Function Parity test Tolerance
sim_es ANN, AAA (with seasonal initial) 1e-10
sim_ssarima AR(1), IMA(1,1) 1e-9
sim_sma order 3 1e-10
sim_oes odds-ratio, inverse-odds-ratio, direct, general 1e-10
sim_ces none, full 1e-9
sim_gum local-level, 2-component 1e-10
OMG.simulate combined-probability 1e-9

Known divergences in OM.simulate / simulate.om

Three R-parity tests on the simulate-from-fit path are currently marked xfail with explanations:

  1. test_om_simulate_latent_matches_r — OM fits on the same y converge to slightly different smoothing parameters between Python (NLopt) and R for small intermittent series. With matched errors via the plug-in trick the latent series differ by the ratio alpha_python / alpha_R. This is the underlying OM-fit-parity issue, not a simulate-path divergence.

  2. test_om_simulate_seasonal_matches_r — multiplicative- seasonal OM (MNM, lags=[1, 12]) produces all-NaN latent on the Python side because the multiplicative seasonal state can collapse to zero under matched errors. R handles the same matrix prep without NaN. Tracked as a Python-side ADAM.simulate multiplicative-seasonal numerics follow-up.

  3. test_om_simulate_arima_matches_r — R's simulate.om on an OM with ARMA orders crashes inside the C++ kernel (element-wise multiplication: incompatible matrix dimensions: 2x1 and 1x1). An upstream R-side bug in simulateADAMCore's matrix prep for OM+ARIMA. Python's OM.simulate handles the same case cleanly (the smoke test passes); byte-equivalence is impossible until the R bug is fixed.

The OMG combined-probability case (test_omg_simulate_combined_probability_matches_r) does pass at atol=1e-9 because the OMG combination formula (omg_link_function) operates on the latent series and the optimizer divergence on a single OM fit gets averaged out when combining two independently-fit sub-models.

Practical guidance

When you want to reproduce a specific simulation across both languages:

  • Pre-generate the error vector in either language (it doesn't matter which — both have access to the other's outputs via CSV / JSON / R_data round-trip).
  • Pass it via the randomizer callable form on both sides.
  • The C++ kernel guarantees byte-equivalence (modulo BLAS rounding in the few code paths that touch lstsq).

When you want reproducible-within-language simulations:

  • Use the language's native seed argument. The output is deterministic given the same seed within that language but doesn't translate across.

Current parity-test tolerances

The Python repo's r_parity tests (test_om_summary_r_comparison.py, test_omg_r_comparison.py, test_adam_summary_r_comparison.py, test_fi_r_comparison.py, test_multicov_r_parity.py, test_coefbootstrap_r_parity.py, test_sim_es_r_parity.py, test_sim_ssarima_r_parity.py, test_sim_sma_r_parity.py, test_sim_oes_r_parity.py, test_sim_ces_r_parity.py, test_sim_gum_r_parity.py, test_om_simulate_r_parity.py) carry the r_parity marker and are deselected in CI by default.

Quantity rtol atol
Coefficients (all scenarios) machine precision —
Log-likelihood, loss, fitted, residuals machine precision —
OM / OMG / backcasting-ADAM vcov / SE / CI 1e-4 1e-6
ADAM optimal / two-stage vcov, SE, CI bounds 2e-2 1e-3
multicov(type="analytical") 5e-2 1e-3
multicov(type="empirical") 5e-2 1e-3
coefbootstrap (stochastic — distributional only) structural
sim_* (plug-in-numbers trick) — 1e-9 to 1e-10
OMG.simulate (plug-in-numbers trick) — 1e-9
OM.simulate (plug-in-numbers trick) xfail (3 cases — see Simulation section above)

Clone this wiki locally