pith. sign in

arxiv: 2604.09375 · v1 · submitted 2026-04-10 · 🧮 math.OC · math.ST· stat.ML· stat.TH

Data-Efficient Non-Gaussian Semi-Nonparametric Density Estimation for Nonlinear Dynamical Systems

Aaron R. Liao , Kenshiro Oguri , Michele D. Carpenter This is my paper

Pith reviewed 2026-05-10 16:43 UTC · model grok-4.3

classification 🧮 math.OC math.STstat.MLstat.TH
keywords semi-nonparametric density estimationnon-Gaussian distributionsnonlinear dynamical systemsHermite polynomialsMonte Carlo integrationquantile estimationLorenz systemconvex relaxation
0
0 comments X

The pith

Semi-nonparametric densities using Hermite polynomials estimate non-Gaussian state distributions in nonlinear dynamics with far fewer samples than Monte Carlo.

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

This paper establishes a method to estimate probability density functions for states in nonlinear dynamical systems that exhibit non-Gaussian behavior. It uses semi-nonparametric densities based on Hermite polynomials that remain positive by construction. Monte Carlo sampling approximates the integrals needed for maximum likelihood fitting of the coefficients, and a convex relaxation provides good starting points. Demonstrated on the Lorenz attractor, the approach captures the density shape and quantiles accurately while using significantly less data than direct Monte Carlo sampling, which matters when each simulation run is expensive.

Core claim

The paper claims that by modeling densities with a Hermite polynomial expansion that guarantees positivity, approximating the maximum likelihood objective via Monte Carlo, and initializing coefficients through convex relaxation, one can obtain accurate non-Gaussian density estimates and quantiles for chaotic systems such as the Lorenz equations using substantially fewer forward propagations than would be required by raw Monte Carlo sampling.

What carries the argument

Semi-nonparametric (SNP) densities constructed from a probabilists' Hermite polynomial basis, with Monte Carlo approximation of the expectation integrals arising in maximum likelihood estimation of the coefficients and a convex relaxation to obtain initial estimates.

Load-bearing premise

Monte Carlo approximations of the integrals in the maximum likelihood estimation and the convex relaxation for initial coefficients work reliably when each evaluation of the nonlinear dynamics is expensive.

What would settle it

Simulate the Lorenz system with a limited number of samples using both the proposed method and direct Monte Carlo, then compare the estimated quantiles against a high-fidelity reference; the claim fails if the SNP method does not produce more accurate quantiles than raw Monte Carlo at the same small sample size.

Figures

Figures reproduced from arXiv: 2604.09375 by Aaron R. Liao, Kenshiro Oguri, Michele D. Carpenter.

Figure 2
Figure 2. Figure 2: Objective values comparison across different polynomial orders [PITH_FULL_IMAGE:figures/full_fig_p006_2.png] view at source ↗
Figure 1
Figure 1. Figure 1: Monte Carlo Point Cloud To obtain a density estimate, first, Ns MC points are generated from the initial distribution. These points are propagated through the Lorenz dynamics and used to solve the convex relaxed SNP optimization. The coefficients from the relaxed problem are then used as initial guesses to the nonlinear SNP optimization problem [PITH_FULL_IMAGE:figures/full_fig_p006_1.png] view at source ↗
Figure 4
Figure 4. Figure 4: SNP density generated from 1000 MC points with [PITH_FULL_IMAGE:figures/full_fig_p007_4.png] view at source ↗
Figure 6
Figure 6. Figure 6: Quantile evaluation comparisons As expected, the 10 trials of 106 samples provide the most accurate estimate, predicting an average of 6.67247% of the distribution lies within the defined box with a very tight spread as shown with the respective box plot. Meanwhile, the SNP density estimates perform well in this scenario, getting very close to the 6.67247% estimate from the 106 MC samples with far fewer sa… view at source ↗
Figure 5
Figure 5. Figure 5: Density estimate with Monte Carlo cloud and constraint box [PITH_FULL_IMAGE:figures/full_fig_p007_5.png] view at source ↗
read the original abstract

Accurate representation of non-Gaussian distributions of quantities of interest in nonlinear dynamical systems is critical for estimation, control, and decision-making, but can be challenging when forward propagations are expensive to carry out. This paper presents an approach for estimating probability density functions of states evolving under nonlinear dynamics using Seminonparametric (SNP), or Gallant-Nychka, densities. SNP densities employ a probabilists' Hermite polynomial basis to model non-Gaussian behavior and are positive everywhere on the support by construction. We use Monte Carlo to approximate the expectation integrals that arise in the maximum likelihood estimation of SNP coefficients, and introduce a convex relaxation to generate effective initial estimates. The method is demonstrated on density and quantile estimation for the chaotic Lorenz system. The results demonstrate that the proposed method can accurately capture non-Gaussian density structure and compute quantiles using significantly fewer samples than raw Monte Carlo sampling.

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, simulated authors' rebuttal, and a circularity audit. Tearing a paper down is the easy half of reading it; the pith above is the substance, this is the friction.

Referee Report

3 major / 2 minor

Summary. The paper proposes a data-efficient method for estimating non-Gaussian probability density functions of states evolving under nonlinear dynamics. It employs semi-nonparametric (SNP) densities using a probabilists' Hermite polynomial basis (positive by construction), approximates the expectation integrals in maximum-likelihood estimation of the coefficients via Monte Carlo sampling, and introduces a convex relaxation to obtain reliable initial estimates. The approach is demonstrated on density and quantile estimation for the chaotic Lorenz system, with the central claim that it accurately captures non-Gaussian structure using significantly fewer forward propagations than direct Monte Carlo sampling.

Significance. If the central claim holds under broader conditions, the work would be significant for uncertainty quantification, estimation, and control in nonlinear systems where each forward propagation is computationally expensive. Strengths include the positivity guarantee of the SNP representation, the use of Monte Carlo for tractable MLE, and the convex relaxation for initialization; these are explicitly credited as enabling data efficiency. The Lorenz demonstration shows promise for capturing tails and modes, but broader validation would be needed to establish impact beyond cheap simulations.

major comments (3)
  1. [Abstract] Abstract and demonstration section: the claim that the method 'accurately capture[s] non-Gaussian density structure and compute[s] quantiles using significantly fewer samples than raw Monte Carlo sampling' is not supported by any quantitative error metrics, convergence rates, L2 or quantile errors, or ablation studies against baselines. Without these, it is impossible to verify whether Monte Carlo approximation error or the convex relaxation introduces bias that affects the central claim.
  2. [Numerical results] Numerical example (Lorenz system): the only test case uses inexpensive forward propagations. The central claim requires that MC-approximated MLE plus convex relaxation remain stable and accurate when each propagation is expensive and the optimizer is limited to a few hundred runs; no such stress test, variance analysis of the MC integrals, or multimodal likelihood check is provided, leaving the data-efficiency claim unverified for the intended regime.
  3. [Method] Method description: the convex relaxation is introduced to generate initial coefficient estimates, yet no analysis or bound is given on how the relaxation error propagates into the final SNP density or quantile estimates, which is load-bearing for the claim of reliable non-Gaussian tail capture.
minor comments (2)
  1. [Method] Notation for the SNP polynomial coefficients and the Monte Carlo sample size should be introduced with explicit symbols and distinguished from the number of forward propagations.
  2. [Numerical results] Figure captions for the Lorenz density plots should include the number of samples used in the proposed method versus raw MC and report the specific quantiles shown.

Simulated Author's Rebuttal

3 responses · 0 unresolved

We thank the referee for the constructive comments, which have helped strengthen the quantitative support and methodological discussion in our work. We have revised the manuscript to incorporate additional error metrics, variance analysis, and numerical sensitivity studies. Our point-by-point responses to the major comments follow.

read point-by-point responses

  1. Referee: [Abstract] Abstract and demonstration section: the claim that the method 'accurately capture[s] non-Gaussian density structure and compute[s] quantiles using significantly fewer samples than raw Monte Carlo sampling' is not supported by any quantitative error metrics, convergence rates, L2 or quantile errors, or ablation studies against baselines. Without these, it is impossible to verify whether Monte Carlo approximation error or the convex relaxation introduces bias that affects the central claim.

    Authors: We agree that explicit quantitative metrics strengthen the central claim. In the revised manuscript, we have added L2 norm errors (computed against a high-sample reference density), quantile errors at levels 0.05/0.25/0.75/0.95, and convergence plots of these errors versus sample size. Ablation comparisons to raw Monte Carlo and to the method without the convex initialization are now included in the numerical results section, along with error bars from repeated runs. These additions confirm that neither the Monte Carlo approximation nor the relaxation introduces appreciable bias at the sample counts used. revision: yes

  2. Referee: [Numerical results] Numerical example (Lorenz system): the only test case uses inexpensive forward propagations. The central claim requires that MC-approximated MLE plus convex relaxation remain stable and accurate when each propagation is expensive and the optimizer is limited to a few hundred runs; no such stress test, variance analysis of the MC integrals, or multimodal likelihood check is provided, leaving the data-efficiency claim unverified for the intended regime.

    Authors: The Lorenz system is a standard benchmark for nonlinear chaotic dynamics that produce non-Gaussian state distributions. Although the forward model is inexpensive to evaluate in simulation, the method itself is formulated to operate with only a few hundred propagations, which directly addresses the expensive-propagation regime. In the revision we have added (i) an empirical variance analysis of the Monte Carlo integrals appearing in the MLE objective, (ii) a discussion of how the convex relaxation supplies reliable starting points that mitigate multimodal likelihood issues, and (iii) a simulated stress test that caps the optimizer at 200 function evaluations while still recovering accurate density and quantile estimates. These elements verify stability under the limited-budget setting targeted by the paper. revision: partial

  3. Referee: [Method] Method description: the convex relaxation is introduced to generate initial coefficient estimates, yet no analysis or bound is given on how the relaxation error propagates into the final SNP density or quantile estimates, which is load-bearing for the claim of reliable non-Gaussian tail capture.

    Authors: We acknowledge that an analytic bound on relaxation-error propagation is not provided. The convex relaxation is used exclusively for initialization and is followed by full non-convex MLE refinement. In the revised manuscript we have included a numerical sensitivity study that perturbs the initial coefficients by amounts comparable to observed relaxation residuals and tracks the resulting change in the final density and quantile estimates; the perturbations are shown to be negligible after optimization. This empirical evidence supports reliable tail capture for the Lorenz example. A general theoretical bound would require additional assumptions on the likelihood surface and is noted as future work, but the practical reliability is now documented. revision: partial

Circularity Check

0 steps flagged

No circularity; derivation is self-contained estimation against external Monte Carlo benchmarks

full rationale

The paper presents SNP density estimation via MC-approximated MLE integrals plus convex relaxation for coefficient initialization, with performance claims evaluated by direct comparison to independent raw Monte Carlo sampling on the Lorenz system. No equation or claim reduces a reported result (e.g., quantile accuracy or sample efficiency) to a quantity fitted from the same data by construction. No self-citation is load-bearing for the central method or uniqueness. The approach is presented as a standard estimation procedure whose validity rests on external numerical verification rather than tautological redefinition of inputs.

Axiom & Free-Parameter Ledger

1 free parameters · 2 axioms · 0 invented entities

The central claim rests on the mathematical guarantee that Hermite-polynomial expansions around a Gaussian yield valid positive densities, plus the assumption that Monte Carlo sampling adequately approximates the likelihood integrals for the chosen sample budget.

free parameters (1)
  • SNP polynomial coefficients
    Estimated via maximum likelihood; their number (polynomial order) is a modeling choice that must be selected or cross-validated.
axioms (2)
  • standard math Probabilists' Hermite polynomials multiplied by a Gaussian weight produce a valid probability density that is positive everywhere on the real line
    Invoked by the choice of SNP/Gallant-Nychka representation; this is a standard property of the basis.
  • domain assumption Monte Carlo sampling provides a sufficiently accurate approximation to the expectation integrals arising in the SNP likelihood
    Central to the data-efficient claim; error in this approximation directly affects coefficient quality.

pith-pipeline@v0.9.0 · 5461 in / 1544 out tokens · 59204 ms · 2026-05-10T16:43:34.628035+00:00 · methodology

discussion (0)

Sign in with ORCID, Apple, or X to comment. Anyone can read and Pith papers without signing in.

Reference graph

Works this paper leans on

25 extracted references · 25 canonical work pages

  1. [1]

    An Introduction to Sequential Monte Carlo Methods,

    A. Doucet, N. de Freitas, and N. Gordon, “An Introduction to Sequential Monte Carlo Methods,” inSequential Monte Carlo Methods in Practice, A. Doucet, N. de Freitas, and N. Gordon, Eds. New York, NY: Springer, 2001, pp. 3–14

    work page 2001

  2. [2]

    M. J. Kochenderfer,Decision Making Under Uncertainty: Theory and Application. MIT Press, Jul. 2015

    work page 2015

  3. [3]

    B. W. Silverman,Density Estimation for Statistics and Data Analysis. London: Chapman and Hall, 1986

    work page 1986

  4. [4]

    Multivariate density estimation,

    D. W. Scott, “Multivariate density estimation,” 2015

    work page 2015

  5. [5]

    Remarks on Some Nonparametric Estimates of a Density Function,

    M. Rosenblatt, “Remarks on Some Nonparametric Estimates of a Density Function,”The Annals of Mathematical Statistics, vol. 27, no. 3, pp. 832–837, Sep. 1956

    work page 1956

  6. [6]

    On Estimation of a Probability Density Function and Mode,

    E. Parzen, “On Estimation of a Probability Density Function and Mode,”The Annals of Mathematical Statistics, vol. 33, no. 3, pp. 1065–1076, Sep. 1962

    work page 1962

  7. [7]

    A Higher-Order Autonomous Navigation Filter For Nonlinear Dynamics And Non-Gaussian Dis- tributions,

    A. R. Liao, K. Oguri, and M. Carpenter, “A Higher-Order Autonomous Navigation Filter For Nonlinear Dynamics And Non-Gaussian Dis- tributions,” inAAS Guidance, Navigation, and Control conference, Breckenridge, CO, Feb. 2026

    work page 2026

  8. [8]

    Non-Gaussian Dis- tribution Steering in Nonlinear Dynamics with Conjugate Unscented Transformation,

    D. C. Qi, K. Oguri, P. Singla, and M. R. Akella, “Non-Gaussian Dis- tribution Steering in Nonlinear Dynamics with Conjugate Unscented Transformation,” Oct. 2025

    work page 2025

  9. [9]

    Chance-Constrained Gaussian Mixture Steering to a Terminal Gaussian Distribution,

    N. Kumagai and K. Oguri, “Chance-Constrained Gaussian Mixture Steering to a Terminal Gaussian Distribution,” in2024 IEEE 63rd Conference on Decision and Control (CDC), Dec. 2024, pp. 2207– 2212

    work page 2024

  10. [10]

    On the Representation of Statistical Frequency by a Series,

    F. Y . Edgeworth, “On the Representation of Statistical Frequency by a Series,”Journal of the Royal Statistical Society, vol. 70, no. 1, pp. 102–106, 1907

    work page 1907

  11. [11]

    Uber die Darstellung willkurlicher Functione,

    C. V . L. Charlier, “Uber die Darstellung willkurlicher Functione,” Arkiv for Matematik, Astronomi och Fysik, vol. 2, pp. 1–35, 1905

    work page 1905

  12. [12]

    Gram–Charlier densities,

    E. Jondeau and M. Rockinger, “Gram–Charlier densities,”Journal of Economic Dynamics and Control, vol. 25, no. 10, pp. 1457–1483, Oct. 2001

    work page 2001

  13. [13]

    The Conditions Under Which Gram- Charlier and Edgeworth Curves are Positive Definite and Unimodal,

    D. E. Barton and K. E. Dennis, “The Conditions Under Which Gram- Charlier and Edgeworth Curves are Positive Definite and Unimodal,” Biometrika, vol. 39, no. 3/4, pp. 425–427, 1952

    work page 1952

  14. [14]

    Multidirectional Gaussian Mixture Mod- els for Nonlinear Uncertainty Propagation,

    V . Vittaldev and R. Russell, “Multidirectional Gaussian Mixture Mod- els for Nonlinear Uncertainty Propagation,”Computer Modeling in Engineering & Sciences, vol. 111, no. 1, pp. 83–117, 2016

    work page 2016

  15. [15]

    Nonlinear Propagation of Orbit Uncertainty Using Non-Intrusive Polynomial Chaos,

    B. A. Jones, A. Doostan, and G. H. Born, “Nonlinear Propagation of Orbit Uncertainty Using Non-Intrusive Polynomial Chaos,”Journal of Guidance, Control, and Dynamics, vol. 36, no. 2, pp. 430–444, 2013

    work page 2013

  16. [16]

    Nonlinear Bayesian estimation using Gaussian sum approximations,

    D. Alspach and H. Sorenson, “Nonlinear Bayesian estimation using Gaussian sum approximations,”IEEE Transactions on Automatic Control, vol. 17, no. 4, pp. 439–448, Aug. 1972

    work page 1972

  17. [17]

    Xiu,Numerical methods for stochastic computations: a spectral method approach

    D. Xiu,Numerical methods for stochastic computations: a spectral method approach. Princeton university press, 2010

    work page 2010

  18. [18]

    Risken,The F okker-Planck Equation: Methods of Solution and Applications, ser

    H. Risken,The F okker-Planck Equation: Methods of Solution and Applications, ser. Springer Series in Synergetics, H. Haken, Ed. Berlin, Heidelberg: Springer, 1996, vol. 18

    work page 1996

  19. [19]

    Semi-Nonparametric Maximum Likelihood Estimation,

    A. R. Gallant and D. W. Nychka, “Semi-Nonparametric Maximum Likelihood Estimation,”Econometrica, vol. 55, no. 2, pp. 363–390, 1987

    work page 1987

  20. [20]

    Seminonparametric Estimation of Conditionally Constrained Heterogeneous Processes: Asset Pricing Applications,

    A. R. Gallant and G. Tauchen, “Seminonparametric Estimation of Conditionally Constrained Heterogeneous Processes: Asset Pricing Applications,”Econometrica, vol. 57, no. 5, pp. 1091–1120, 1989

    work page 1989

  21. [21]

    Econometric Estimators and the Edgeworth Approxi- mation,

    J. D. Sargan, “Econometric Estimators and the Edgeworth Approxi- mation,”Econometrica, vol. 44, no. 3, pp. 421–448, 1976

    work page 1976

  22. [22]

    Multivariate moments expansion den- sity: Application of the dynamic equicorrelation model,

    T.-M. ˜N´ıguez and J. Perote, “Multivariate moments expansion den- sity: Application of the dynamic equicorrelation model,”Journal of Banking & Finance, vol. 72, pp. S216–S232, Nov. 2016

    work page 2016

  23. [23]

    C. P. Robert and G. Casella,Monte Carlo Statistical Methods, ser. Springer Texts in Statistics. New York, NY: Springer, 2004

    work page 2004

  24. [24]

    Ascent-based Monte Carlo expectation– maximization,

    B. S. Caffo, W. Jank, and G. L. Jones, “Ascent-based Monte Carlo expectation– maximization,”Journal of the Royal Statistical Society: Series B (Statistical Methodology), vol. 67, no. 2, pp. 235–251, 2005

    work page 2005

  25. [25]

    A. B. Owen,Monte Carlo theory, methods and examples. Stanford, 2013

    work page 2013