Weak Moment Methods for Statistical Inference
with an Application to Robust Estimation
April 26, 2026
A companion paper [1] develops a generalised probabilistic framework in which a probability law is represented by a tempered distribution \(T\in\mathcal{S}'(\mathbb{R})\)—on the same footing as a density or a characteristic function—and information is extracted by pairing \(T\) with a positive Schwartz kernel \(\varphi\in\mathcal{S}(\mathbb{R})\) that acts as a measurement instrument rather than as part of the law; it shows that the resulting weak moments of all orders exist unconditionally. The present paper turns this structure into a concrete methodology for statistical inference.
We develop estimation strategies based on weak moments, weak characteristic functions, weak cumulants, and regularised reconstruction of the underlying density via Tikhonov inversion of the multiplication operator \(M_\varphi\colon f\mapsto f\varphi\). A key feature of this programme is that parametric inference proceeds directly from weak expectations—moments, transforms, or cumulants—without requiring reconstruction of the underlying density; reconstruction is an additional route, useful when density-level inference is the goal. Each strategy extends a classical procedure to heavy-tailed models in which the corresponding classical object does not exist.
The central contribution is to show that weak moment estimators are automatically locally robust in the sense of Hampel. The score function of a weak moment estimator has the form \(x^j\varphi(x)\) minus the corresponding theoretical weak moment; it is bounded and redescending for every positive Schwartz kernel, its influence function admits a closed form, and its gross error sensitivity is finite in every identifiable parametric model—without ad hoc truncation. The kernel thus plays the role of Huber’s tuning constant, but instead of a post-hoc modification of the estimator it forms part of the chosen characterisation of the model—the prescription of how the law is observed. The arbitrariness inherent in any such tuning is thereby not removed but relocated to a more model-related—and, we conjecture, in some cases more natural—part of the inference.
The programme is worked out in detail for the Cauchy location model (where no classical moment estimator exists), for a Student \(t_3\) location–scale model, for a bivariate Cauchy location model illustrating the multivariate extension, and for a bivariate \(t_3\) location–scale model that reveals the breakdown of the MLE scale estimate under contamination. In each case we compute the influence function, the gross error sensitivity, and the asymptotic variance in closed form, and compare the resulting estimators numerically with the median, the MLE, the Huber \(M\)-estimator, and the Tukey biweight, under both the correctly specified model and under contamination. Although the present paper focuses on parametric models, the reconstruction route is inherently non-parametric and opens a path to weak density estimation without parametric assumptions; we outline this direction in Section sec:sec:sec:discussion?.
Moment-based methods are among the oldest tools in statistics, yet they fail for heavy-tailed models—Cauchy, low-degree Student \(t\), stable laws—whose classical moments do not exist. In a companion paper [1] we developed a framework in which a probability law is carried by a tempered distribution \(T\) and probed by a Schwartz kernel \(\varphi\)—the kernel an instrument for extracting information, not part of the law—with expectations defined via \(\mathbb{E}_{T,\varphi}[\psi]=\langle T,\psi\varphi\rangle\). Because \(\varphi\in\mathcal{S}(\mathbb{R})\) ensures \(x^j\varphi\in\mathcal{S}(\mathbb{R})\), weak moments \({}^{(\varphi)}m_j=\langle T,x^j\varphi\rangle\) exist for every \(j\) and every pair. Labouriau [1] develops the algebra (additivity of weak cumulants, a weak CLT) and the analysis (a hierarchy of uniqueness theorems for the weak moment problem). It closes with a minimal illustration: even the first weak moment yields a consistent estimator for the Cauchy location parameter.
The present paper turns that illustration into a full methodology. Its two objectives are: (i) to develop estimation procedures—weak moment matching, transform-based and cumulant-based methods, regularised reconstruction of the density—inside the weak framework; and (ii) to show that these estimators are automatically locally robust: their influence function is bounded, their gross error sensitivity finite, and their score redescending, all inherited from the decay of the kernel \(\varphi\) with no ad hoc truncation. A key feature is that parametric inference proceeds directly from weak expectations without requiring reconstruction of the underlying density; reconstruction is an additional route, useful when density-level inference is the goal. Although the present paper focuses on parametric models, the reconstruction pathway is inherently non-parametric and opens a direction toward weak density estimation without parametric assumptions (see Section 7). The central methodological point is that inference in the weak framework is fundamentally direct: parameters are identified and estimated from weak expectations without reconstructing the underlying density; reconstruction is only required when the density itself is the inferential target—a secondary, self-contained route in this paper, separable from its direct-inference core and a natural candidate for separate development.
Beyond densities and likelihoods. Two features of the distributional representation deserve emphasis at the outset. First, it does not presuppose a density: because the law is carried by a tempered distribution \(T\), the construction applies unchanged to models that possess no density, and hence no likelihood—non-dominated families, such as a moving atom superimposed on a continuous background, for which likelihood-based inference has no starting point. Second, the representation is structurally stable: by the classical structure theorem every tempered distribution is a finite sum of derivatives of ordinary (polynomially bounded, continuous) functions (Strichartz [2], §6.3; cf. [1]), so singular laws—point masses, jumps—are differentiated regular functions rather than pathologies, and the kernel returns them to bounded, smooth weak quantities. Since the tempered distribution \(T\) characterises the law, this structural stability is a property of the model itself—of the probability measures in play—and not of the estimation method.
Why robustness is the natural home of weak moment methods. At first sight, weak moment estimators might appear to be merely a device for handling heavy-tailed models. We argue that the connection to robust statistics is deeper. First, in the classical programme of Hampel and Huber [3], [4], one starts from an efficient but fragile score (typically the MLE score) and truncates or redescends it to bound the influence function. The resulting estimator depends on a tuning constant. In the weak framework the score \(x^j\varphi(x)-{}^{(\varphi)}m_j(\theta)\) is automatically bounded and redescending, and the “tuning” is the kernel \(\varphi\)—part of the chosen characterisation of the model, not a post-hoc modification. Second, finite gross error sensitivity holds for essentially any positive Schwartz kernel and any parametric family, regardless of the tail behaviour of \(f_\theta\). In particular, the framework produces locally robust estimators for the Cauchy distribution, where the classical Hampel programme has no efficient starting point. Third, the weak framework supplies a family of estimators parametrised by the moment order \(j\), the kernel \(\varphi\), and GMM weighting. Optimising over this family is the analogue of the classical Hampel-optimal estimation problem.
Notation. Following [1], weak objects carry a superscript \(\varphi\): \({}^{(\varphi)}m_j\) for weak moments, \({}^{(\varphi)}\phi(t)\) for the weak characteristic function, \({}^{(\varphi)}K(t)\) for the weak cumulant generating function, \({}^{(\varphi)}\kappa_j\) for weak cumulants. Classical objects carry no superscript. The superscript is a label, not a power.
Organisation. Section 2 recalls the minimal setup from [1] and introduces parametric weak models with their empirical counterparts. Section 3 develops the estimation strategies: direct methods (weak moment matching, transform-based and cumulant methods) and density reconstruction via regularised inversion. Section 4—the centrepiece—places weak estimators inside the theory of robust statistics. Section 5 works out the Cauchy and Student \(t_3\) examples in detail and illustrates the multivariate extension through elliptically contoured models. Section 6 presents four Monte Carlo studies—univariate Cauchy location, univariate \(t_3\) location–scale, bivariate Cauchy location, and bivariate \(t_3\) location–scale—comparing weak moment estimators with classical and robust benchmarks under both the clean model and contamination. Section 7 concludes.
We assume the reader is familiar with the framework of [1]; this section fixes notation and records, without proofs, the properties we shall use.
A distribution–kernel pair is \((T,\varphi)\) with \(T\in\mathcal{S}'(\mathbb{R})\) and \(\varphi\in\mathcal{S}(\mathbb{R})\); the weak expectation of a polynomially bounded \(\psi\) with \(\psi\varphi\in\mathcal{S}(\mathbb{R})\) is \(\mathbb{E}_{T,\varphi}[\psi]:=\langle T,\psi\varphi\rangle\). Throughout, \(T\) represents the probability law and \(\varphi\) is the fixed kernel through which it is measured; following [1] the kernel is an instrument, not part of the law, so every weak quantity below is a measurement of the single underlying law \(T\). In the case that a probability density exists and is sufficiently regular (\(T=T_f\), \(f\ge 0\)), \(\mathbb{E}_{T_f,\varphi}[\psi]=\int\psi(x)\varphi(x)f(x)\,dx\). The weak moments \({}^{(\varphi)}m_j:=\mathbb{E}_{T,\varphi}[x^j]\), the weak characteristic function \({}^{(\varphi)}\phi(t):=\mathbb{E}_{T,\varphi}[e^{itx}]\), and the weak cumulants \({}^{(\varphi)}\kappa_j:=(1/i^j)(d^j/dt^j)\log{}^{(\varphi)}\phi(t)|_{t=0}\) are all well defined for every pair and every \(j\). The uniqueness theorems of [1] guarantee that the weak moment sequence determines the underlying distribution under mild conditions on \(\varphi\) (Gaussian kernels: via Hermite completeness; positive Schwartz kernels: via a Carleman condition; exponential-decay kernels: via Denjoy–Carleman quasi-analyticity). In the density case, weak expectations also determine the kernel-weighted distribution function and the regularised density \(g=f\varphi\) (see [1]); this fact underpins the reconstruction strategy of Section 3.3.
Definition 1 (Generalised random variable). A generalised random variable is a random object whose distribution is specified by a distribution–kernel pair \((T,\varphi)\), and whose expectations are defined by \(\mathbb{E}_{T,\varphi}[\psi]=\langle T,\psi\varphi\rangle\). See [1] for a full treatment.
Remark 1 (Interpretation in the density case). When \((T,\varphi)\) arises from a density \(f\), a generalised random variable can be understood as an ordinary random variable \(X\sim f\) observed through the kernel \(\varphi\): the observations \(X_1,\ldots,X_n\) are in practice ordinary draws from \(f\), and the kernel enters only through the inferential functionals \(\psi(X_i)\varphi(X_i)\). That is, the “generalised” character resides in the expectation operator, not in the sampling mechanism.
Thus, [1] establishes the probabilistic and analytic foundations of the weak framework, while the present paper develops its statistical and inferential consequences.
Let \(\Theta\subseteq\mathbb{R}^p\) be open, let \(\varphi\in\mathcal{S}(\mathbb{R})\) with \(\varphi>0\), and let \(\{T_\theta\}_{\theta\in\Theta}\) be a family of tempered distributions such that each \((T_\theta,\varphi)\) is a probability pair. We call \(\{(T_\theta,\varphi)\}_{\theta\in\Theta}\) a parametric weak model. In the density case the weak moments are \({}^{(\varphi)}m_j(\theta)=\int x^j\varphi(x)f_\theta(x)\,dx\).
The kernel is an instrument for extracting information about the model from the tempered distribution \(T\): different kernels yield different weak moments and hence different estimators. Structural identifiability—injectivity of \(\theta\mapsto T_\theta\)—is a property of the distributional family and is preserved by any positive kernel under the uniqueness results of [1]. In estimation, however, one works with a finite collection of weak moments or transforms, and local identifiability is expressed through the full-rank condition on the Jacobian \(G(\theta)\) of Proposition 4. The kernel therefore affects not only efficiency but also the operative identifiability of the estimating equations.
Remark 2 (Three levels of identifiability). Three notions of identifiability are relevant in this framework: (i) structural identifiability of \(\{T_\theta\}\), a property of the distributional family alone; (ii) identifiability from the full weak moment sequence, which under positive kernels inherits the uniqueness theorems of [1]; and (iii) local identifiability from the chosen estimating equations, governed by \(\operatorname{rank}G(\theta)\). Level (i) is model-intrinsic; level (ii) connects to the weak moment problem; level (iii) is procedure-dependent and involves the kernel through the moment map \(\theta\mapsto{}^{(\varphi)}m_j(\theta)\).
Let \(X_1,\ldots,X_n\) be i.i.d.generalised random variables with distribution \((T_\theta,\varphi)\) in the sense of Definition 1. The empirical weak expectation is \[\label{eq:emp-weak} \widehat{\mathbb{E}}_{}^{(\varphi)}n[\psi] :=\frac{1}{n}\sum_{i=1}^n\psi(X_i)\,\varphi(X_i).\tag{1}\] The factor \(\varphi(X_i)\) inside the average is essential: it matches the theoretical weak expectation \(\int\psi\varphi f_\theta\,dx\), ensures that \(\psi\varphi\) is bounded (hence in \(L^2(f_\theta)\)), and produces automatic downweighting of observations far from the kernel’s centre. By the law of large numbers, \(\widehat{\mathbb{E}}_{}^{(\varphi)}n[\psi]\to\mathbb{E}_{T_\theta,\varphi}[\psi]\) a.s., and the CLT gives asymptotic normality with variance \(\mathrm{Var}_{f_\theta}(\psi\varphi) =\mathbb{E}_{T_\theta,\varphi^2}[\psi^2]-(\mathbb{E}_{T_\theta,\varphi}[\psi])^2\), which is finite unconditionally.
The general estimation principle is therefore: match empirical and theoretical weak expectations, \(\widehat{\mathbb{E}}_{}^{(\varphi)}n[\psi]\approx\mathbb{E}_{T_\theta,\varphi}[\psi]\), over a family of test functions \(\psi\). Specialising to \(\psi_j(x)=x^j\) gives weak moment estimators; to \(\psi_t(x)=e^{itx}\) gives transform-based estimators; to derivatives of \(\log\phi_{}^{(\varphi)}\theta\) gives weak cumulant methods.
It is useful to read the kernel as a level of observational resolution. The weak expectation \(\mathbb{E}_{T,\varphi}[\psi]=\langle T,\psi\varphi\rangle\) is a measurement of the law \(T\) at the resolution set by \(\varphi\): the kernel localises the probe \(\psi\) to the region where \(\varphi\) is appreciable and damps its tails, so that even unbounded or non-integrable features of \(T\) are seen through a finite aperture. Two limiting regimes make this reading precise.
Stability. Tempered distributions act continuously on Schwartz space, so the measurement is stable under smooth changes of the instrument: if \(\varphi_j\to\varphi\) in \(\mathcal{S}(\mathbb{R})\) then \(\langle T,\psi\varphi_j\rangle\to\langle T,\psi\varphi\rangle\) ([1]). Weak moments, weak characteristic functions, and the estimators built from them therefore depend continuously on the observational kernel; the dependence on \(\varphi\) is not a formal artefact but a controlled, intrinsically regularised one.
Concentration. If the kernel is replaced by an approximate identity \(\rho_{\varepsilon}(x)=\varepsilon^{-1}\rho((x-x_0)/\varepsilon)\) with \(\rho\in\mathcal{S}(\mathbb{R})\) and \(\rho_\varepsilon\to\delta_{x_0}\), the measurement concentrates at \(x_0\); in the density case \(\mathbb{E}_{T_f,\rho_\varepsilon}[\psi]\to\psi(x_0)f(x_0)\) at points of continuity, recovering pointwise values of the weighted density.
Recovery of the classical object. Sending the kernel to the constant \(1\) removes the aperture. For the Gaussian kernel \(\varphi_\sigma(x)=e^{-x^2/(2\sigma^2)}\) this is the limit \(\sigma\to\infty\), and whenever the classical expectation \(\int\psi\,dF\) exists one recovers it, \(\mathbb{E}_{T_f,\varphi_\sigma}[\psi]\to\int\psi\,dF\), by dominated convergence; the general \(\varphi\to1\) statement is [1]. Inferentially, the weak estimators of Section 3 are therefore classical estimators observed at finite resolution: as \(\sigma\to\infty\) they reduce to their classical counterparts—the sample mean, the empirical characteristic function—wherever those exist, and they remain well defined (bounded score, redescending influence, finite variance) precisely in the heavy-tailed and singular cases where the classical objects break down. The bandwidth \(\sigma\) thus interpolates between classical efficiency (large \(\sigma\)) and robustness (small \(\sigma\)), a trade-off taken up in Section 4.
Remark 3 (Structure and a microlocal perspective). The resolution reading meshes with the fine structure of tempered distributions. By the structure theorem, \(T\) is a finite sum \(\sum_\alpha D^\alpha f_\alpha\) of derivatives of polynomially bounded continuous functions (Strichartz [2], §6.3; Hörmander [5]; see also [1]); pairing with \(\psi\varphi\) transfers the derivatives onto the instrument, \[\langle T,\psi\varphi\rangle =\sum_\alpha(-1)^{|\alpha|}\!\int f_\alpha\,D^\alpha(\psi\varphi),\] so the finite aperture \(\varphi\) and its derivatives are exactly what render the singular (differentiated) part of the law as finite scalar measurements. A finer question—in which directions the singularities of \(T\) lie—belongs to the microlocal analysis of distributions (singular support, wave front set), where a kernel acts as a device that attenuates selected microlocal features; this geometric direction is developed in the companion work on transversality [6] and is not pursued here.
The estimation strategies developed in this section fall into two groups. The first three—weak moment matching (Section 3.1), transform-based methods, and cumulant methods (Section 3.2)—estimate parameters directly from weak data, without reconstructing the underlying density. The fourth—regularised reconstruction (Section 3.3)—recovers the density itself and is needed only when density-level inference is the goal.
Fix a set of moment orders \(\mathcal{J}\subseteq\mathbb{N}\) and a positive definite weighting matrix \(W\). The weak moment estimator is \[\label{eq:wm-def} \hat{\theta}_n(\mathcal{J},W) :=\arg\min_{\theta\in\Theta}\; g_n(\theta)^\top W\,g_n(\theta),\tag{2}\] where \(g_n(\theta):=(\hat{}^{(\varphi)}m_j-{}^{(\varphi)}m_j(\theta))_{j\in\mathcal{J}}\) and \(\hat{}^{(\varphi)}m_j:=n^{-1}\sum_i X_i^j\varphi(X_i)\). When \(|\mathcal{J}|=p\) the system is just-identified and the estimator solves \(\hat{}^{(\varphi)}m_j={}^{(\varphi)}m_j(\hat{\theta})\); when \(|\mathcal{J}|>p\) a two-step GMM with \(W=\hat{S}^{-1}\) is asymptotically optimal. The system is well-posed for every parametric weak model because weak moments exist unconditionally.
Proposition 4 (Asymptotics). Under standard regularity—\(m(\cdot;\varphi)\in C^1\), the Jacobian \(G(\theta):=\partial {}^{(\varphi)}m(\theta)/\partial\theta\) of full rank, \(\theta_0\) interior to \(\Theta\)—the weak moment estimator is consistent and \[\sqrt{n}(\hat{\theta}_n-\theta_0) \xrightarrow{d} \mathcal{N}\bigl(0,\,V(\theta_0;W)\bigr),\] where \(V=(G^\top WG)^{-1}G^\top WS\,WG(G^\top WG)^{-1}\), \(S_{jk}(\theta)=\mathbb{E}_{T_\theta,\varphi^2}[x^{j+k}] -{}^{(\varphi)}m_j(\theta)\,{}^{(\varphi)}m_k(\theta)\), and the optimal weight is \(W=S(\theta_0)^{-1}\).
Proof. The result follows from classical GMM/\(M\)-estimation theory once the regularity conditions are verified; the key observation is that the moment functions \(x^j\varphi(x)\) lie in \(L^2(f_\theta)\) for every Schwartz kernel, which ensures all required moment and smoothness conditions. A self-contained verification is given in Appendix 8.
Remark 5 (Connection with the distributional CLT). The asymptotic normality established in Proposition 4 can also be understood through the distributional CLT of [1]. In the density case, the kernel-weighted function \(h=\varphi f\) is itself a probability density with all moments finite (since \(\varphi\in\mathcal{S}(\mathbb{R})\) and \(f\in L^1\)), so the classical CLT applies to \(h\)-distributed observations. The weak moment \({}^{(\varphi)}m_j(\theta)\) is precisely the \(j\)-th moment of \(h\), and the GMM asymptotics of Proposition 4 can be viewed as a specialisation of this classical CLT to the kernel-weighted measure. This provides a second, independent justification of the asymptotic normality of weak moment estimators.
The empirical weak characteristic function \(\hat{}^{(\varphi)}\phi_n(t):=n^{-1}\sum_i e^{itX_i}\varphi(X_i)\) consistently estimates \(\phi_{}^{(\varphi)}\theta(t)\), since \(e^{itx}\varphi(x)\) is bounded. Estimation can proceed by minimising a discrepancy \(\int|\hat{}^{(\varphi)}\phi_n(t)-\phi_{}^{(\varphi)}\theta(t)|^2\,w(t)\,dt\) for a weight \(w\in L^1\). Weak cumulants provide a third route: matching empirical and theoretical weak cumulants, obtained from the moment–cumulant recursion applied to \(\{{}^{(\varphi)}m_j\}\). Both approaches share the fundamental property that all ingredients are finite and smooth, even when the classical counterparts are not.
The preceding methods estimate parameters directly from weak data, without reconstructing the density. We now describe a complementary route for settings where the density itself is the inferential target. It is logically self-contained—separable from the direct-inference core of the paper, and may be read or skipped on its own—and its full development (operator inversion, regularisation rates, minimax theory) is extensive enough to warrant separate treatment; we confine ourselves here to the elements needed to show that the weak framework supports density-level inference in principle.
The weak characteristic function \({}^{(\varphi)}\phi_{T_f,\varphi}(t)\) is the Fourier transform of \(g:=f\varphi=M_\varphi f\), where \(M_\varphi\colon L^2(\mathbb{R})\to L^2(\mathbb{R})\) is pointwise multiplication by \(\varphi\). This reconstruction strategy builds on results established in [1]: weak expectations determine the kernel-weighted distribution function ([1]), and in the density case the regularised density \(g=f\varphi\) can be stably recovered by Tikhonov inversion of the multiplication operator \(M_\varphi\) ([1]). We now develop this into an estimation procedure and refine the convergence analysis.
Proposition 6. \(M_\varphi\) is bounded, self-adjoint, and injective. Its inverse is densely defined but unbounded.
Proof. Boundedness: \(\|M_\varphi f\|_2\le\|\varphi\|_\infty\|f\|_2\). Self-adjointness: \(\varphi\) is real. Injectivity: \(\varphi>0\) and \(\varphi f=0\) a.e.imply \(f=0\). Since \(\inf_x\varphi(x)=0\), \(M_\varphi^{-1}\) is unbounded. ◻
Naive inversion \(f=g/\varphi\) amplifies noise where \(\varphi\) is small. The Tikhonov-regularised reconstruction \[\label{eq:tikh} (R_\lambda g)(x) =\frac{\varphi(x)}{\varphi(x)^2+\lambda}\,g(x), \qquad\lambda>0,\tag{3}\] is the unique minimiser of \(\|M_\varphi h-g\|_{L^2}^2+\lambda\|h\|_{L^2}^2\).
Remark 7 (Three layers of reconstruction). The reconstruction pathway has three layers: (i) weak data determine kernel-weighted distributional information [1]; (ii) in the density case, this reduces to the inverse problem \(g=M_\varphi f\), whose regularised solution is analysed in [1]; and (iii) the present theorem (Theorem 8 below) refines the convergence with a source-condition rate and frames the reconstruction as an estimation strategy. Layers (i)–(ii) are foundational; layer (iii) is inferential.
Theorem 8. If \(f\in L^2(\mathbb{R})\) and \(g=f\varphi\), then \(\|R_\lambda g-f\|_{L^2}\to 0\) as \(\lambda\to 0^+\). Under the source condition \(f=\varphi^\nu h\) with \(\nu\in(0,2)\), the rate is \(\|R_\lambda g-f\|_{L^2}\le c_\nu\lambda^{\nu/2}\|h\|_{L^2}\).
Proof. The \(L^2\)-consistency as \(\lambda\downarrow 0\) is established in [1]; the source-condition rate is the new contribution of the present theorem. Under \(f=\varphi^\nu h\), substituting \(s=\varphi(x)/\sqrt\lambda\) shows \(|R_\lambda g-f|^2\le\lambda^\nu c_\nu^2|h|^2\), and integration concludes. ◻
The kernel mediates a fundamental trade-off: rapid decay of \(\varphi\) ensures existence of weak moments for a broad class of distributions (including heavy-tailed ones) but makes reconstruction more ill-posed; slow decay facilitates reconstruction but regularises less. The analyst can choose the kernel—a degree of freedom not usually available in inverse problems, where the forward operator is dictated by physics.
Remark 9 (Non-parametric character of reconstruction). The reconstruction machinery of this section is stated in a parametric context for consistency with the rest of the paper, but it is inherently non-parametric: no parametric form of \(f\) is assumed in the Tikhonov inversion 3 . The procedure takes empirical weak data, recovers \(g=\varphi f\), and inverts \(M_\varphi\) by regularisation, without parametric constraints. A systematic non-parametric treatment—including minimax rates under smoothness or source conditions on \(f\), and the trade-off between tail regularisation and ill-posedness as a function of \(\varphi\)—is deferred to future work (see Section 7).
This section is the centrepiece of the paper. We show that every weak moment estimator is a locally robust \(M\)-estimator in the sense of Hampel, with bounded influence function, finite gross error sensitivity, and a redescending score—all inherited from the kernel.
Let \(T\colon\mathcal{F}\to\mathbb{R}^p\) be a statistical functional. The influence function at \(F\) is \(\mathrm{IF}(x;T,F):=\lim_{\varepsilon\downarrow 0} \varepsilon^{-1}[T((1-\varepsilon)F+\varepsilon\delta_x)-T(F)]\), and the gross error sensitivity is \(\gamma^{*}(T,F):=\sup_x\|\mathrm{IF}(x;T,F)\|\). An estimator is \(B\)-robust at \(F\) if \(\gamma^{*}<\infty\). For an \(M\)-estimator with score \(\psi\), \[\label{eq:IF-M} \mathrm{IF}(x;T_\psi,F_\theta) =-[M(\theta)]^{-1}\psi(x;\theta), \qquad M(\theta):=\int\frac{\partial\psi}{\partial\theta}\,dF_\theta,\tag{4}\] and the asymptotic variance is \(V=M^{-1}QM^{-\top}\) with \(Q=\int\psi\psi^\top\,dF_\theta\). Within this class, \(B\)-robustness reduces to boundedness of \(\psi\).
Fix \(\varphi\in\mathcal{S}(\mathbb{R})\) with \(\varphi>0\) and \(j\in\mathbb{N}\). Define the score \[\label{eq:score} \psi_j(x;\theta;\varphi) :=x^j\varphi(x)-{}^{(\varphi)}m_j(\theta).\tag{5}\] The weak moment estimating equation \(n^{-1}\sum_i\psi_j(X_i;\hat{\theta};\varphi)=0\) is the \(M\)-estimation equation with this score.
Proposition 10 (Bounded score). For every \(\theta\), \(j\), and \(\varphi\in\mathcal{S}(\mathbb{R})\), the score \(\psi_j(\cdot;\theta;\varphi)\) is bounded and rapidly decreasing.
Proof. \(x^j\varphi(x)\in\mathcal{S}(\mathbb{R})\) because \(\mathcal{S}\) is closed under polynomial multiplication. ◻
Proposition 11 (Influence function). For a scalar parameter \(\theta\) and a single moment order \(j\) (the just-identified case), if \(\partial_\theta {}^{(\varphi)}m_j(\theta)\ne 0\), the influence function of the weak moment estimator is \[\label{eq:IF-wm} \mathrm{IF}(x;T_{\psi_j},F_\theta) =\frac{{}^{(\varphi)}m_j(\theta)-x^j\varphi(x)}{\partial_\theta {}^{(\varphi)}m_j(\theta)}.\qquad{(1)}\]
Proof. The score depends on \(\theta\) only through the constant \({}^{(\varphi)}m_j(\theta)\), so \(M(\theta)=-\partial_\theta {}^{(\varphi)}m_j(\theta)\). Applying 4 gives the result. ◻
When \(\theta\in\mathbb{R}^p\), or when several moment orders are stacked, the scalar derivative \(\partial_\theta{}^{(\varphi)}m_j\) is replaced by the Jacobian \(G(\theta)=\partial{}^{(\varphi)}m(\theta)/\partial\theta\), and the influence function takes the matrix form of Proposition 12; equation ?? is its scalar, just-identified specialisation.
Corollary 1 (Gross error sensitivity). \(\gamma^{*}(T_{\psi_j},F_\theta) =\sup_x|x^j\varphi(x)-{}^{(\varphi)}m_j(\theta)| /|\partial_\theta {}^{(\varphi)}m_j(\theta)|<\infty\).
Corollary 1 is, in our view, the central observation of this paper: every weak moment estimator in a parametric weak model with identifiable parameter has finite gross error sensitivity, with no conditions on the tails of \(f_\theta\), no tuning, and no truncation—that is, it is locally (\(B\)-)robust. This concerns the influence function and gross error sensitivity, not the breakdown point or maximum bias, which we address separately in Section 4.5.
Since \(x^j\varphi(x)\to 0\) as \(|x|\to\infty\), the influence function ?? redescends to a finite constant for large \(|x|\)—the hallmark of redescending \(M\)-estimators in the sense of [7].
For a Gaussian kernel \(\varphi_\sigma(x)=e^{-x^2/(2\sigma^2)}\), the bandwidth \(\sigma\) controls the redescent scale: small \(\sigma\) gives strong downweighting (like a Tukey biweight with a small tuning constant); large \(\sigma\) gives slow redescent and recovers the classical estimator in the limit \(\sigma\to\infty\). The kernel thus plays the same role as the tuning constant in a classical redescending \(M\)-estimator, but arises from the model specification rather than being introduced solely for robustness.
As with all redescending \(M\)-estimators, the estimating equation may have multiple roots. In practice this is handled by initialising the root-finder at a robust preliminary estimator (e.g.the sample median), by choosing \(\sigma\) large enough to cover the plausible parameter range, or by using a GMM formulation with several moment orders. The connection between the kernel’s decay, the identifiability region, and the multiplicity of roots deserves emphasis: kernels with faster decay reduce the gross error sensitivity but may also narrow the effective identifiability region, since the score \(\psi_j(x;\theta;\varphi)\) becomes flatter over a wider range of \(\theta\). This trade-off between robustness and operative identifiability is governed by the kernel bandwidth and is closely related to the three levels of identifiability discussed in Remark 2.
Let \(\mathcal{J}=\{j_1,\ldots,j_K\}\) and define the vector score \(\Psi(x;\theta;\varphi):=(x^{j_k}\varphi(x) -{}^{(\varphi)}m_{j_k}(\theta))_{k=1}^K\).
Proposition 12 (GMM influence function). The influence function of the GMM estimator with weighting \(W\) is \[\label{eq:IF-GMM} \mathrm{IF}(x;\hat{\theta}_n,F_\theta) =(G^\top WG)^{-1}G^\top W\,\Psi(x;\theta;\varphi),\qquad{(2)}\] where \(G=-\partial {}^{(\varphi)}m(\theta)/\partial\theta\). It is bounded, and \(\gamma^{*}<\infty\) whenever \(G\) has full rank.
Proof. Standard GMM influence function formula [7] applied to the bounded score \(\Psi\). ◻
The optimal \(W=S(\theta)^{-1}\) minimises the asymptotic variance but not necessarily the GES; optimising \(W\) subject to a GES constraint is the weak analogue of Hampel’s optimality problem.
Local robustness (bounded IF, finite GES) is automatic in the weak framework. Global robustness (breakdown point, maximum bias) depends more delicately on the kernel and on the model. For a Gaussian kernel, observations with \(|X_i|\gg\sigma\) contribute essentially zero to the empirical weak moment, so very extreme contamination has negligible effect; but intermediate contamination within the kernel’s effective support can affect the estimator. A systematic study of the maximum bias curve and finite-sample breakdown as functions of \(\sigma\) and \(\mathcal{J}\) is left for future work; the simulation study in Section 6 provides empirical evidence of good performance under moderate contamination.
We collect the principal conclusions:
Every weak moment estimator has a bounded, rapidly decreasing score (Proposition 10).
Its influence function has the closed form ?? , bounded whenever the parameter is identifiable.
Its gross error sensitivity is finite in any parametric weak model (Corollary 1), without conditions on the tails of \(f_\theta\).
The estimator is redescending, with redescent rate governed by the kernel.
The asymptotic variance has a closed form in terms of weak moments (Proposition 4) and is finite unconditionally.
GMM with multiple weak moments retains all properties and typically improves efficiency.
The overall picture: the weak framework provides, by construction, a family of locally robust \(M\)-estimators whose design parameter is a kernel rather than a truncation constant.
Consider the Cauchy location model \(f(x;\mu)=[\pi(1+(x-\mu)^2)]^{-1}\), where no classical moment exists. We use the Gaussian kernel \(\varphi_\sigma(x)=e^{-x^2/(2\sigma^2)}\) and the first weak moment \(j=1\).
Using the Fourier representation \(\hat{f}(\cdot;\mu)(t)=e^{i\mu t-|t|}\) and Parseval’s theorem, \({}^{(\varphi)}m_0(\mu;\sigma)=\int\varphi_\sigma f(\cdot;\mu)\,dx\) reduces to a Voigt-profile integral expressible via the Faddeeva function at \((\mu+i)/(\sigma\sqrt{2})\). The first weak moment satisfies \({}^{(\varphi)}m_1(\mu;\sigma)=\mu\,{}^{(\varphi)}m_0(\mu;\sigma) +\mathcal{R}(\mu;\sigma)\) with a correction \(\mathcal{R}\) computable from the same special function. In practice all integrals are evaluated by adaptive quadrature.
The normalised first weak moment \({}^{(\varphi)}m_1/{}^{(\varphi)}m_0\) is strictly increasing in \(\mu\) for \(|\mu|<\mu^*(\sigma)\); for \(\sigma=3\), \(\mu^*\approx 8\), which covers any practical range.
The score \(\psi(x;\mu;\sigma)=x\varphi_\sigma(x)-{}^{(\varphi)}m_1(\mu;\sigma)\) is a bounded, smooth, redescending function of \(x\). By Proposition 11, \[\mathrm{IF}(x;\hat{\mu},f(\cdot;\mu)) =\frac{{}^{(\varphi)}m_1(\mu;\sigma)-x\varphi_\sigma(x)}{\partial_\mu {}^{(\varphi)}m_1(\mu;\sigma)}.\] At \(\mu=0\) (by symmetry, \({}^{(\varphi)}m_1(0)=0\)) and \(\sup_x x\,e^{-x^2/(2\sigma^2)}=\sigma/\sqrt{e}\), \[\gamma^{*}=\frac{\sigma/\sqrt{e}}{|\partial_\mu {}^{(\varphi)}m_1(0;\sigma)|}.\] For \(\sigma=3\), \(\gamma^{*}\approx 5.5\)—larger than the median’s GES (\(\pi/2\approx 1.57\)) but comparable to a Huber estimator with moderate tuning. The asymptotic variance is \(V(\mu;\sigma)=S(\mu;\sigma)/(\partial_\mu {}^{(\varphi)}m_1)^2\) with \(S=\mathbb{E}_{T_\mu,\varphi^2}[x^2]-({}^{(\varphi)}m_1)^2\); at \(\sigma=3\), \(\mu=0\), \(V\approx 3.1\) versus the median’s \(\pi^2/4\approx 2.47\) (relative efficiency \(\approx 0.80\)). GMM with multiple moments closes most of this gap.
Consider the Student \(t_\nu\) location–scale model with \(\nu=3\) (classical mean and variance exist; third and higher moments diverge). We use \(\varphi_\sigma\) with \(\sigma=3\) and two weak moments \(j\in\{1,2\}\) for the two parameters \((\mu,s)\).
By Proposition 12, the influence function is \[\mathrm{IF}(x;(\hat{\mu},\hat{s}),f(\cdot;\mu,s)) =(G^\top WG)^{-1}G^\top W \begin{pmatrix} x\varphi_\sigma(x)-{}^{(\varphi)}m_1(\mu,s;\sigma)\\ x^2\varphi_\sigma(x)-{}^{(\varphi)}m_2(\mu,s;\sigma) \end{pmatrix},\] where \(G=-\partial({}^{(\varphi)}m_1,{}^{(\varphi)}m_2)/\partial(\mu,s)\). It is bounded whenever \(G\) is invertible (generic). The GES is finite, and the redescent rate is governed by \(\sigma\) exactly as in the Cauchy case.
We briefly illustrate the extension of weak moment methods to \(\mathbb{R}^d\) through elliptically contoured models. Let \(X\in\mathbb{R}^d\) have density \[f(x;\mu,\Sigma) =|\Sigma|^{-1/2}\,g\!\bigl((x-\mu)^\top\Sigma^{-1}(x-\mu)\bigr),\] where \(\mu\in\mathbb{R}^d\), \(\Sigma\) is positive definite, and \(g\) is a radial profile. Classical moments may fail to exist when \(g\) is heavy-tailed. Let \(\varphi_\sigma(x)=\exp(-\|x\|^2/(2\sigma^2))\) be an isotropic Gaussian kernel. Then for each coordinate \(k=1,\ldots,d\), the weak first moment \[{}^{(\varphi)}m_{1,k}(\mu,\Sigma;\sigma) =\int_{\mathbb{R}^d}x_k\,\varphi_\sigma(x)\,f(x;\mu,\Sigma)\,dx\] is well defined for all parameter values, regardless of the tail behaviour of \(f\).
A natural estimator of the location vector \(\mu\) is obtained by matching empirical and theoretical weak moments: \[\widehat\mu_n\;\text{ solves }\; \frac{1}{n}\sum_{i=1}^n X_i\,\varphi_\sigma(X_i) ={}^{(\varphi)}m_1(\mu,\Sigma;\sigma),\] where \({}^{(\varphi)}m_1\) denotes the vector of weak first moments. When \(\Sigma\) is known, this yields a \(d\)-dimensional estimating equation; when \(\Sigma\) is unknown, additional weak second moments can be included in a GMM formulation.
The corresponding score function \[\psi(x;\mu,\Sigma;\sigma) =x\,\varphi_\sigma(x)-{}^{(\varphi)}m_1(\mu,\Sigma;\sigma)\] is bounded and rapidly decreasing in \(\|x\|\), so the influence function is bounded and redescending componentwise. In particular, the local robustness properties established in Section 4 extend directly to the multivariate setting.
As a concrete instance, consider the multivariate Cauchy with density \[f(x;\mu) =\frac{\Gamma\!\bigl(\tfrac{d+1}{2}\bigr)}{\pi^{(d+1)/2}} \bigl(1+\|x-\mu\|^2\bigr)^{-(d+1)/2},\] for which no classical moment of any order exists. With the isotropic Gaussian kernel, the weak first-moment vector satisfies \[{}^{(\varphi)}m_{1,k}(\mu;\sigma) =\int_{\mathbb{R}^d}x_k\,e^{-\|x\|^2/(2\sigma^2)} \frac{\Gamma\!\bigl(\tfrac{d+1}{2}\bigr)}{\pi^{(d+1)/2}} \bigl(1+\|x-\mu\|^2\bigr)^{-(d+1)/2}\,dx.\] By the translation structure \(f(x;\mu)=f_0(x-\mu)\), a change of variables gives \({}^{(\varphi)}m_{1,k}(\mu;\sigma) =\mu_k\,{}^{(\varphi)}m_0(\mu;\sigma)+\mathcal{R}_k(\mu;\sigma)\), where \({}^{(\varphi)}m_0\) is the weak zeroth moment (a convolution of a Gaussian kernel with a heavy-tailed radial profile—the multivariate analogue of the classical Voigt integral that arises when Gaussian and Lorentzian line shapes are convolved [8]) and \(\mathcal{R}_k\) is a correction computable by adaptive quadrature. In practice, specialising to \(d=2\) or \(d=3\) and evaluating the integrals numerically is straightforward. The uniqueness of the weak moment sequence for Gaussian kernels on \(\mathbb{R}^d\) is guaranteed by [1], so the estimation programme is well-posed.
This example illustrates that weak moment methods extend naturally to \(\mathbb{R}^d\): the kernel ensures existence of moments and boundedness of scores, while the estimation structure remains identical to the univariate case. The main additional challenge is computational, as the number of moments required for joint estimation of \((\mu,\Sigma)\) grows with the dimension. A Monte Carlo illustration for \(d=2\) is given in Section 6.3.
Remark 13 (Elliptical laws without explicit densities). The density formulation above is used for concreteness. The weak framework also applies to elliptically contoured laws specified through their characteristic functions or distributional representations, including cases where a density is unavailable in closed form. For instance, symmetric stable laws on \(\mathbb{R}^d\) with characteristic function \(\exp(-\|t\|^\alpha)\), \(\alpha\in(0,2)\), have no known closed-form density for general \(d\) and \(\alpha\), yet the weak expectation is still defined through the distributional pairing, and the same empirical weak estimating equations can be used whenever the corresponding theoretical weak moments or transforms are computable.
The preceding examples are heavy-tailed but dominated. We now turn to a model that admits no dominating measure—and hence no likelihood—to show that weak moment estimation applies without change. Consider the location family \[P_\theta=\tfrac12\,\delta_\theta+\tfrac12\,N(0,1), \qquad\theta\in\mathbb{R},\] a point mass of weight \(\tfrac12\) at the unknown location \(\theta\) superimposed on a standard Gaussian background. Since \(P_\theta(\{\theta\})=\tfrac12\) for every \(\theta\), a common dominating measure would need an atom at every point of \(\mathbb{R}\), which no \(\sigma\)-finite measure possesses; the family is therefore non-dominated. There is no common density, no likelihood, and the Fisher–Rao machinery does not apply. As a tempered distribution, however, \(T_\theta=\tfrac12\delta_\theta+\tfrac12 N(0,1)\in\mathcal{S}'(\mathbb{R})\) is perfectly regular, and its weak moments against any Schwartz kernel are finite; the Gaussian kernel determines \(T_\theta\) uniquely (the determinacy example of [1]).
With the Gaussian kernel \(\varphi_\sigma(x)=e^{-x^2/(2\sigma^2)}\) the weak first moment is available in closed form; the symmetric Gaussian background contributes nothing to it, so \[{}^{(\varphi)}m_1(\theta;\sigma) =\tfrac12\,\theta\,e^{-\theta^2/(2\sigma^2)}, \qquad \partial_\theta{}^{(\varphi)}m_1(\theta;\sigma) =\tfrac12\,e^{-\theta^2/(2\sigma^2)} \Bigl(1-\tfrac{\theta^2}{\sigma^2}\Bigr).\] The map \(\theta\mapsto{}^{(\varphi)}m_1\) is strictly increasing on \((-\sigma,\sigma)\), so the atom location is identified from the first weak moment whenever \(|\theta|<\sigma\), and \(\widehat\theta\) solving \(n^{-1}\sum_i X_i\varphi_\sigma(X_i)={}^{(\varphi)}m_1(\widehat\theta;\sigma)\) is \(\sqrt n\)-consistent there. Additional weak moments in a GMM (Section 4.4) enlarge the identifiable range beyond \(|\theta|<\sigma\).
The estimator is the \(M\)-estimator with the bounded score \(\psi(x;\theta;\sigma)=x\varphi_\sigma(x)-{}^{(\varphi)}m_1(\theta;\sigma)\) of Section 4, so its closed forms apply verbatim. By Proposition 11, \[\mathrm{IF}(x;\widehat\theta,P_\theta) =\frac{{}^{(\varphi)}m_1(\theta;\sigma)-x\varphi_\sigma(x)}{\tfrac12 e^{-\theta^2/(2\sigma^2)}(1-\theta^2/\sigma^2)},\] bounded and redescending in \(x\) because \(\sup_x|x\varphi_\sigma(x)| =\sigma/\sqrt e\). The gross error sensitivity and asymptotic variance are \[\gamma^{*}(\theta;\sigma) =\frac{\sigma/\sqrt e+|{}^{(\varphi)}m_1(\theta;\sigma)|}{|\partial_\theta{}^{(\varphi)}m_1(\theta;\sigma)|}, \qquad V(\theta;\sigma) =\frac{\mathbb{E}_{T_\theta,\varphi^2}[x^2]-{}^{(\varphi)}m_1(\theta;\sigma)^2}{\bigl(\partial_\theta{}^{(\varphi)}m_1(\theta;\sigma)\bigr)^2},\] with \(\mathbb{E}_{T_\theta,\varphi^2}[x^2] =\tfrac12\theta^2e^{-\theta^2/\sigma^2} +\tfrac12\sigma^3/(2+\sigma^2)^{3/2}\). Both are finite on \(|\theta|<\sigma\) and grow as \(|\theta|\uparrow\sigma\), where the first weak moment ceases to identify \(\theta\); Table 1 gives representative values for \(\sigma=3\).
| \(\theta\) | \({}^{(\varphi)}m_1\) | \(\partial_\theta{}^{(\varphi)}m_1\) | \(\GES\) | \(V(\theta;\sigma)\) |
|---|---|---|---|---|
| \(0.0\) | \(0.000\) | \(0.500\) | \(3.64\) | \(1.48\) |
| \(0.5\) | \(0.247\) | \(0.479\) | \(4.31\) | \(1.88\) |
| \(1.0\) | \(0.473\) | \(0.420\) | \(5.45\) | \(3.36\) |
| \(1.5\) | \(0.662\) | \(0.331\) | \(7.50\) | \(7.38\) |
| \(2.0\) | \(0.801\) | \(0.222\) | \(11.78\) | \(20.44\) |
The point is structural: although the model has no density and no likelihood, the weak first moment supplies a smooth estimating equation with a bounded, redescending influence function and a closed-form asymptotic variance—the same inferential apparatus used for the dominated examples, applied unchanged to a non-dominated family.
We present Monte Carlo comparisons (\(2{,}000\) replications) of weak moment estimators against classical benchmarks and robust estimators, under both the correctly specified model and under contamination.
True model: Cauchy\((\mu,1)\) with \(\mu=2\), Gaussian kernel \(\sigma=3\). Estimators: WM (single moment \(j=1\)), GMM-I (identity weighting, \(\mathcal{J}=\{1,2\}\)), GMM-2S (two-step optimal, ridge \(\lambda=0.10\)), Median, MLE (Newton from median), Huber (\(k=1.345\)), Tukey biweight (\(c=4.685\)). Huber and Tukey constants are calibrated for \(95\%\) Gaussian asymptotic efficiency.
| WM | GMM-I | GMM-2S | Median | MLE | Huber | Tukey | ||||||||
| \(n\) | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE |
| 50 | 0.01 | 0.29 | \(-\)0.01 | 0.30 | 0.02 | 0.28 | 0.00 | 0.23 | 0.00 | 0.18 | 0.01 | 0.33 | 0.01 | 0.31 |
| 100 | 0.01 | 0.20 | 0.00 | 0.20 | 0.00 | 0.18 | 0.00 | 0.16 | 0.00 | 0.13 | 0.00 | 0.23 | 0.00 | 0.22 |
| 500 | 0.00 | 0.09 | 0.00 | 0.09 | 0.00 | 0.08 | 0.00 | 0.07 | 0.00 | 0.06 | 0.00 | 0.10 | 0.00 | 0.10 |
| 1000 | 0.00 | 0.06 | 0.00 | 0.06 | 0.00 | 0.05 | 0.00 | 0.05 | 0.00 | 0.04 | 0.00 | 0.07 | 0.00 | 0.07 |
| 5000 | 0.00 | 0.03 | 0.00 | 0.03 | 0.00 | 0.02 | 0.00 | 0.02 | 0.00 | 0.02 | 0.00 | 0.03 | 0.00 | 0.03 |
Contaminated model: \((1-\varepsilon)\,\text{Cauchy}(\mu,1) +\varepsilon\,\text{Cauchy}(\mu+\delta,1)\) with \(\varepsilon=0.10\), \(\delta=5\).
| WM | GMM-I | GMM-2S | Median | MLE | Huber | Tukey | ||||||||
| \(n\) | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE | Bias | RMSE |
| 50 | 0.08 | 0.31 | 0.16 | 0.34 | 0.12 | 0.33 | 0.16 | 0.30 | 0.18 | 0.27 | 0.14 | 0.34 | 0.09 | 0.31 |
| 100 | 0.09 | 0.24 | 0.19 | 0.29 | 0.06 | 0.22 | 0.16 | 0.25 | 0.18 | 0.22 | 0.13 | 0.26 | 0.07 | 0.23 |
| 500 | 0.07 | 0.12 | 0.18 | 0.20 | 0.05 | 0.10 | 0.16 | 0.18 | 0.18 | 0.19 | 0.12 | 0.14 | 0.06 | 0.12 |
| 1000 | 0.07 | 0.10 | 0.18 | 0.19 | 0.04 | 0.08 | 0.15 | 0.16 | 0.18 | 0.19 | 0.12 | 0.13 | 0.05 | 0.09 |
| 5000 | 0.07 | 0.07 | 0.18 | 0.18 | 0.04 | 0.05 | 0.15 | 0.16 | 0.18 | 0.18 | 0.12 | 0.12 | 0.05 | 0.06 |
The contamination experiment reverses the ranking: the MLE is now the worst performer, the median is biased, and the GMM-2S estimator achieves the smallest bias and RMSE at \(n\ge 100\), comparable to the Tukey biweight—without any hand-tuning beyond the kernel bandwidth.
True model: \(t_3(\mu,s)\) with \((\mu,s)=(0,1)\), Gaussian kernel \(\sigma=3\), \(\mathcal{J}=\{1,2\}\). Comparisons: WM-GMM-2S, MLE, Mean/SD, Median/MAD, Tukey biweight (\(c=4.685\)). We report RMSE under the clean model and under scale contamination \(0.9\,t_3(0,1)+0.1\,t_3(0,5)\).
| Clean \(t_3\) | Contaminated | ||||||
| Estimator | Param. | \(n=100\) | \(n=500\) | \(n=1000\) | \(n=100\) | \(n=500\) | \(n=1000\) |
| WM-GMM-2S | \(\mu\) | 0.18 | 0.08 | 0.06 | 0.20 | 0.09 | 0.06 |
| \(s\) | 0.19 | 0.09 | 0.06 | 0.22 | 0.10 | 0.07 | |
| MLE | \(\mu\) | 0.16 | 0.07 | 0.05 | 0.19 | 0.09 | 0.06 |
| \(s\) | 0.15 | 0.07 | 0.05 | 0.18 | 0.09 | 0.07 | |
| Mean/SD | \(\mu\) | 0.19 | 0.08 | 0.06 | 0.44 | 0.20 | 0.14 |
| \(s\) | 0.46 | 0.20 | 0.14 | 0.92 | 0.42 | 0.30 | |
| Median/MAD | \(\mu\) | 0.20 | 0.09 | 0.06 | 0.22 | 0.10 | 0.07 |
| \(s\) | 0.24 | 0.11 | 0.08 | 0.30 | 0.14 | 0.10 | |
| Tukey | \(\mu\) | 0.18 | 0.08 | 0.06 | 0.20 | 0.09 | 0.06 |
| \(s\) | 0.20 | 0.09 | 0.06 | 0.24 | 0.11 | 0.08 | |
Under the clean model, the MLE is most efficient; WM-GMM-2S, Tukey, and Median/MAD are within \(10\)–\(20\%\). Under scale contamination, the Mean/SD pair is severely damaged; the weak moment estimator, Tukey biweight, and Median/MAD retain performance, with WM-GMM-2S competitive throughout.
We illustrate the multivariate extension of Section 5.3 in dimension \(d=2\). Data are generated from the bivariate Cauchy location model \[f(x;\mu) =\frac{\Gamma(3/2)}{\pi^{3/2}} \bigl(1+\|x-\mu\|^2\bigr)^{-3/2}, \qquad x\in\mathbb{R}^2,\] with \(\mu=(1,1)^\top\). Here we use the isotropic Gaussian kernel \(\varphi_\sigma(x)=\exp\{-\|x\|^2/(2\sigma^2)\}\) with \(\sigma=3\) and estimate \(\mu\) by matching the vector weak first moment. The estimating equation is solved numerically using the coordinatewise median as starting value. We compare the weak moment (WM) estimator with the coordinatewise median, the spatial median [9], and the multivariate Cauchy MLE.
| WM | Spatial Med. | Coord.Med. | MLE | |||||
| \(n\) | \(\|\text{Bias}\|\) | RMSE | \(\|\text{Bias}\|\) | RMSE | \(\|\text{Bias}\|\) | RMSE | \(\|\text{Bias}\|\) | RMSE |
| Clean model | ||||||||
| 50 | 0.03 | 0.37 | 0.00 | 0.29 | 0.01 | 0.32 | 0.00 | 0.27 |
| 100 | 0.01 | 0.26 | 0.01 | 0.21 | 0.01 | 0.23 | 0.01 | 0.19 |
| 500 | 0.01 | 0.11 | 0.00 | 0.09 | 0.00 | 0.10 | 0.00 | 0.08 |
| 1000 | 0.01 | 0.08 | 0.00 | 0.06 | 0.00 | 0.07 | 0.00 | 0.06 |
| Contaminated: \(0.9\,\mathrm{Cauchy}_2(\mu,I)+0.1\,\mathrm{Cauchy}_2(\mu+\delta,I)\), \(\delta=(5,5)^\top\) | ||||||||
| 50 | 0.10 | 0.37 | 0.21 | 0.39 | 0.23 | 0.44 | 0.07 | 0.30 |
| 100 | 0.11 | 0.26 | 0.20 | 0.31 | 0.22 | 0.34 | 0.07 | 0.21 |
| 500 | 0.11 | 0.15 | 0.20 | 0.22 | 0.22 | 0.25 | 0.07 | 0.11 |
| 1000 | 0.12 | 0.14 | 0.19 | 0.21 | 0.22 | 0.23 | 0.07 | 0.09 |
Under the clean model, the MLE is most efficient; the WM estimator and both medians are within \(20\)–\(40\%\). Under contamination, the ranking changes: the WM estimator outperforms both the spatial and coordinatewise medians in both bias and RMSE. The Cauchy MLE retains good performance because the Cauchy score \(2(x-\mu)/(1+\|x-\mu\|^2)\) is itself bounded and redescending—a well-known property of the Cauchy family [7]—so that the MLE is automatically robust in this particular model. This is consistent with the message of Section 4: for the Cauchy, both the MLE and the WM estimator are redescending \(M\)-estimators; they differ in the source of redescent (likelihood vs kernel).
The bivariate Cauchy example above shows that both the MLE and the WM estimator are robust, owing to the redescending character of the Cauchy score. To exhibit a setting where the MLE does break down, we turn to the bivariate Student \(t_3\) location–scale model: \[f(x;\mu,s) = \frac{\Gamma(5/2)}{\Gamma(3/2)\,(3\pi s^2)} \Bigl(1+\frac{\|x-\mu\|^2}{3s^2}\Bigr)^{-5/2}, \qquad x\in\mathbb{R}^2,\] with \(\mu=(\mu_1,\mu_2)^\top\) and common scale \(s>0\). The parameter vector is \(\theta=(\mu_1,\mu_2,s)\). We use the isotropic Gaussian kernel \(\varphi_\sigma(x)=\exp\{-\|x\|^2/(2\sigma^2)\}\) with \(\sigma=3\), and match three weak moments: the two first-order moments \({}^{(\varphi)}m_{(1,0)}(\theta)\) and \({}^{(\varphi)}m_{(0,1)}(\theta)\) and a sum-of-squares moment \({}^{(\varphi)}m_{(2)}(\theta)=\mathbb{E}_\theta[(X_1^2+X_2^2)\,\varphi_\sigma(X)]\), using a GMM-2S procedure.
The contaminated model is a scale contamination: \(0.9\cdot t_3(\mu,I)+0.1\cdot t_3(\mu,5I)\), which inflates the dispersion of \(10\%\) of the observations without shifting the centre. The location MLE should remain reasonable, since the \(t_3\) score for \(\mu\) is bounded (though not strongly redescending); but the scale MLE, governed by the score \(\partial\log f/\partial s\), is sensitive to large observations and should be severely biased upward.
We compare four strategies: the WM-GMM-2S estimator; the \(t_3\) MLE (iterative reweighting); the sample mean and standard deviation (Mean/SD); and the coordinatewise median with MAD-based scale (Med/MAD). Each estimator returns an estimate of \((\mu,s)\); Table 6 reports the RMSE of \(\hat{\mu}\) (Euclidean norm of the bias vector) and of \(\hat{s}\) separately, over \(2{,}000\) replications.
| WM-GMM-2S | MLE | Mean/SD | Med/MAD | ||
| \(n\) | Param | RMSE | RMSE | RMSE | RMSE |
| Clean model | |||||
| 100 | \(\mu\) | 0.18 | 0.17 | 0.25 | 0.20 |
| \(s\) | 0.09 | 0.08 | 0.81 | 0.17 | |
| 500 | \(\mu\) | 0.08 | 0.08 | 0.11 | 0.09 |
| \(s\) | 0.04 | 0.03 | 0.75 | 0.14 | |
| 1000 | \(\mu\) | 0.06 | 0.05 | 0.08 | 0.06 |
| \(s\) | 0.03 | 0.03 | 0.74 | 0.14 | |
| Contaminated: \(0.9\,t_3(\mu,I)+0.1\,t_3(\mu,5I)\) | |||||
| 100 | \(\mu\) | 0.17 | 0.18 | 0.46 | 0.21 |
| \(s\) | 0.10 | 0.19 | 2.23 | 0.29 | |
| 500 | \(\mu\) | 0.08 | 0.08 | 0.20 | 0.09 |
| \(s\) | 0.06 | 0.16 | 2.19 | 0.27 | |
| 1000 | \(\mu\) | 0.06 | 0.06 | 0.14 | 0.07 |
| \(s\) | 0.06 | 0.16 | 2.18 | 0.26 | |
Under the clean model, the MLE and WM-GMM-2S estimators are comparable for both location and scale; the Mean/SD pair is heavily biased in scale because the sample standard deviation estimates the marginal standard deviation \(s\sqrt{\nu/(\nu-2)}=s\sqrt{3}\approx 1.73\) rather than the scale parameter \(s=1\), producing a systematic bias of approximately \(0.73\). Under scale contamination, the critical finding is in the scale parameter: the MLE scale RMSE remains at approximately \(0.16\) even at \(n=1{,}000\), indicating that the MLE scale estimate does not converge to the true scale under contamination. In contrast, the WM-GMM-2S scale RMSE drops to \(0.06\) at \(n=1{,}000\), converging at the parametric rate. The Med/MAD estimator is robust but less efficient than WM-GMM-2S.
This example complements the bivariate Cauchy study: there, both the MLE and the WM estimator were robust (from different sources); here, the MLE scale estimator breaks down under contamination while the WM estimator retains full performance. The difference is structural: the Cauchy score for location is naturally redescending, but the \(t_3\) score for scale is not—it grows without bound for large observations. The weak moment score, in contrast, inherits redescent from the kernel regardless of the underlying family.
Finally we illustrate the non-dominated model \(P_\theta=\tfrac12\delta_\theta+\tfrac12N(0,1)\) of Section 5.4, with \(\theta=1\) and Gaussian kernel \(\sigma=3\). Because the model has no likelihood there is no MLE; we compare the weak moment estimator (WM, \(j=1\)) with the sample median and with the method-of-moments estimator \(\widehat\theta=2\bar X\) (since \(\mathbb{E}_{P_\theta}[X]=\theta/2\)), over \(2{,}000\) replications, under the clean model and under \(10\%\) contamination by \(N(\theta+5,1)\).
| WM | Median | MoM \(2\bar X\) | ||||
| \(n\) | Bias | RMSE | Bias | RMSE | Bias | RMSE |
| Clean model | ||||||
| \(50\) | \(+0.02\) | \(0.27\) | \(-0.03\) | \(0.11\) | \(+0.01\) | \(0.25\) |
| \(100\) | \(+0.00\) | \(0.18\) | \(-0.01\) | \(0.05\) | \(-0.01\) | \(0.17\) |
| \(500\) | \(+0.00\) | \(0.08\) | \(+0.00\) | \(0.00\) | \(+0.00\) | \(0.08\) |
| \(1000\) | \(-0.00\) | \(0.06\) | \(+0.00\) | \(0.00\) | \(-0.00\) | \(0.05\) |
| \(10\%\) contamination by \(N(\theta+5,1)\) | ||||||
| \(50\) | \(+0.11\) | \(0.28\) | \(-0.01\) | \(0.05\) | \(+1.12\) | \(1.23\) |
| \(100\) | \(+0.10\) | \(0.21\) | \(-0.00\) | \(0.02\) | \(+1.10\) | \(1.16\) |
| \(500\) | \(+0.09\) | \(0.13\) | \(+0.00\) | \(0.00\) | \(+1.10\) | \(1.12\) |
| \(1000\) | \(+0.09\) | \(0.11\) | \(+0.00\) | \(0.00\) | \(+1.10\) | \(1.10\) |
Two points emerge. First, the weak estimator behaves as predicted: its clean-model RMSE matches the closed-form asymptotic standard error of Table 1, and under contamination its bias remains bounded—indeed it redescends, the contribution of an outlier at \(x\) vanishing as \(|x|\to\infty\) because \(x\varphi_\sigma(x)\to0\). Second, the method-of-moments estimator, built on the unbounded statistic \(\bar X\), is overwhelmed by a \(10\%\) contamination. The median is near-exact in this idealised model because the atom carries half the mass and lies at the median; the weak estimator trades that sharpness for smoothness, differentiability, and the closed-form asymptotics shared with the rest of the paper, and it continues to apply when the atom is smeared by measurement noise—so that exact ties, and with them the median’s special behaviour, disappear.
The methodology developed in this paper rests on turning the distributional framework of [1] into a working programme for statistical inference. From the single device of replacing densities by distribution–kernel pairs and defining expectations through the pairing \(\langle T,\psi\varphi\rangle\), we obtain a coherent approach to estimation in heavy-tailed models: parameters can be estimated directly from weak moments, weak characteristic functions, or weak cumulants, without reconstructing the underlying density; and the resulting estimators are automatically locally robust, with bounded influence function, finite gross error sensitivity, and a redescending score—all inherited from the rapid decay of the kernel \(\varphi\).
The central message is that the kernel simultaneously defines the model, regularises the moments, and shapes the influence function. This unification of moment-based and robust inference through generalised probability is, in our view, the principal conceptual contribution. In the classical programme of Hampel and Huber [3], [4], [7], robust estimators require careful tuning of truncation or redescent constants; here the kernel provides this tuning as a structural component of the model, not a post-hoc modification.
The instrument is the tuning. This dual role is clarified by the viewpoint of the companion framework [1], in which \(\varphi\) is not part of the probability law but the instrument through which the law \(T\) is observed. The present paper shows that this same instrument is the robustness tuning: the gross error sensitivity, the redescent scale, and the asymptotic efficiency are all fixed by \(\varphi\) (Corollary 1, Proposition 4). Choosing how to observe the law and choosing how robust the estimator should be are therefore one and the same act—in contrast to the classical programme, where the model is specified first and a tuning constant is grafted on afterwards to bound the influence function.
Where weak estimators are most useful. The Cauchy location example shows that weak moment methods yield consistent, robust estimators where no classical moment-based method exists. In that example, both the MLE and the WM estimator are robust—the former because the Cauchy score is naturally redescending, the latter because of the kernel. The bivariate \(t_3\) location–scale example then reveals the complementary picture: the MLE scale estimator breaks down under contamination (its RMSE does not converge), while the WM-GMM estimator retains full performance at the parametric rate. This demonstrates that automatic local robustness via the kernel is a genuine advantage whenever the likelihood score is not naturally bounded. More generally, the framework is most valuable in settings where tuning is undesirable: unlike the Huber or Tukey estimators [4], [10], whose performance depends on correctly calibrated constants, weak moment estimators depend on a kernel that is part of the model specification.
The role of the distributional CLT. The distributional CLT established in [1] provides a second justification for the asymptotic normality of weak moment estimators (Remark 5). Because the kernel-weighted density \(h=\varphi f\) has all moments finite, the classical CLT applies to \(h\)-distributed observations, and the GMM asymptotics of Proposition 4 can be viewed as a specialisation. This connection underscores that weak inference is not a departure from classical statistics but a regularised extension of it.
Several questions remain open and define a natural research programme.
Optimal kernel selection. Among all positive Schwartz kernels, which minimises the asymptotic variance subject to a gross error sensitivity constraint? This is the weak analogue of Hampel’s optimality problem [7]. For a Gaussian kernel \(\varphi_\sigma\), the bandwidth \(\sigma\) mediates a trade-off between efficiency (large \(\sigma\)) and robustness (small \(\sigma\)); optimising this trade-off in specific parametric families is a natural first step.
Global robustness. The automatic local robustness established in Section 4 concerns the influence function and gross error sensitivity. A systematic study of the maximum bias curve and finite-sample breakdown point for weak moment estimators, as functions of the kernel and the moment set \(\mathcal{J}\), is needed. For Gaussian kernels, the simulation evidence in Section 6 suggests good performance under moderate contamination, but a rigorous analysis would require tools from the theory of maximum bias curves [11].
Efficient GMM and the Cramér–Rao bound. For a fixed family and kernel, what is the best achievable asymptotic variance, and how does it compare with the parametric efficiency bound? The GMM-2S estimator with optimal weighting is efficient within the class of weak moment estimators, but its relationship to the Fisher information remains to be clarified. This is connected to the classical theory of semiparametric efficiency [12] adapted to the weak setting.
Multivariate extensions. The framework extends naturally to \(\mathbb{R}^d\), as illustrated for elliptically contoured models in Section 5.3. In [1], the weak moment problem is shown to have a unique solution for Gaussian kernels on \(\mathbb{R}^d\) via the completeness of the multivariate Hermite basis. This guarantees that the estimation programme of Section 3.1 extends in principle to \(\mathbb{R}^d\): the weak moment sequence determines the distribution, and empirical weak moments provide consistent estimating equations. The main open challenge is computational: the number of moments of order at most \(k\) in \(\mathbb{R}^d\) grows as \(\binom{k+d}{d}\), and optimal selection of the moment set \(\mathcal{J}\) becomes non-trivial. Regularisation strategies analogous to penalised GMM may be needed. Extensions to dependent observations (time series, regression) are likewise natural but unexplored.
Non-parametric density estimation. The reconstruction strategy of Section 3.3 extends naturally to a fully non-parametric setting: observe empirical weak expectations or the empirical weak characteristic function, recover \(g=\varphi f\), and invert \(M_\varphi\) by Tikhonov or spectral regularisation. Sharp minimax rates under smoothness or source conditions on \(f\), and the trade-off between tail regularisation and ill-posedness as a function of \(\varphi\), remain to be established. This programme sits at the intersection of classical kernel smoothing [13], deconvolution problems [14], and robust heavy-tail regularisation. The distinctive feature of the weak approach is that the forward operator \(M_\varphi\) is chosen by the analyst, rather than dictated by an observation model—a degree of freedom that could be exploited to optimise the bias–variance trade-off.
Connections with other frameworks. The automatic robustness of weak moment estimators connects to the broader programme of redescending \(M\)-estimators initiated by Hampel [3] and developed by Beaton and Tukey [10]. The key difference is that in the classical approach, redescent is imposed by modifying the score after the model is specified, whereas here it arises from the model itself through the kernel. There is also a connection to the theory of inverse problems in statistics [13], [15]: the reconstruction route of Section 3.3 is an inverse problem with a chosen forward operator, linking weak inference to regularisation theory. A systematic exploration of these connections is expected to yield further insights into the interplay between distributional modelling, robustness, and inverse methods.
These questions define a coherent programme growing out of the present framework, linking distributional probability, robust statistics, inverse problems, and weighted approximation theory.
We verify the standard regularity conditions for GMM consistency and asymptotic normality (see, e.g., [12], [16], [17]) and show that they are satisfied automatically by the Schwartz property of the kernel.
Setting. Let \(\varphi\in\mathcal{S}(\mathbb{R})\) with \(\varphi>0\), let \(\mathcal{J}=\{j_1,\ldots,j_K\}\) be a set of moment orders, and write \(g(x;\theta):=\bigl(x^{j_k}\varphi(x) -{}^{(\varphi)}m_{j_k}(\theta)\bigr)_{k=1}^K\) for the moment function vector. The GMM estimator minimises \(Q_n(\theta):=g_n(\theta)^\top W\,g_n(\theta)\) with \(g_n(\theta):=n^{-1}\sum_{i=1}^n g(X_i;\theta)\).
Condition 1 (Identification). The population criterion \(Q(\theta):=\mathbb{E}[g(X;\theta)]^\top W\) \(\mathbb{E}[g(X;\theta)]\) is uniquely minimised at \(\theta_0\) whenever the Jacobian \(G(\theta_0):=\partial\mathbb{E}[g(X;\theta)]/\partial\theta \big|_{\theta_0}\) has full column rank and \(W\) is positive definite. This is assumed in the statement of Proposition 4.
Condition 2 (Uniform law of large numbers). For each \(\theta\), \(|g(x;\theta)|\le|x^{j_K}\varphi(x)| +\max_k|{}^{(\varphi)}m_{j_k}(\theta)|\). Since \(x^{j_K}\varphi(x)\in\mathcal{S}(\mathbb{R})\), it is bounded and hence \(\mathbb{E}[|g(X;\theta)|^2]<\infty\) for every \(\theta\). Combined with the continuity of \(\theta\mapsto g(x;\theta)\) and compactness of a neighbourhood of \(\theta_0\), this gives a pointwise (in fact uniform) law of large numbers for \(g_n(\theta)\) by Theorem 2.4 of [17].
Condition 3 (Asymptotic normality of the moment vector). By the multivariate CLT, \(\sqrt{n}\,g_n(\theta_0)\xrightarrow{d} \mathcal{N}(0,S)\), where \(S_{jk}=\mathrm{Cov}\bigl(X^{j_j}\varphi(X),\,X^{j_k}\varphi(X)\bigr)\). Each entry of \(S\) has the form \(\mathbb{E}_{T_\theta,\varphi^2}[x^{j_j+j_k}] -{}^{(\varphi)}m_{j_j}(\theta)\,{}^{(\varphi)}m_{j_k}(\theta)\), which is a weak moment with respect to the pair \((T_\theta,\varphi^2)\) and hence exists unconditionally (since \(\varphi^2\in\mathcal{S}(\mathbb{R})\) whenever \(\varphi\in\mathcal{S}(\mathbb{R})\)).
Condition 4 (Smoothness of the moment map). The assumed regularity \(m(\cdot;\varphi)\in C^1(\Theta)\) ensures that \(G(\theta)\) is continuous. Differentiation under the integral sign is justified because \(|(\partial/\partial\theta_\ell)\,g(x;\theta)|\le |(\partial/\partial\theta_\ell)\,{}^{(\varphi)}m_{j_k}(\theta)|\), which is integrable by the smoothness assumption.
Conclusion. With Conditions 1–4 in hand, the standard GMM theorem [17] yields consistency, and the delta-method argument of [16] yields the asymptotic distribution \(\sqrt{n}(\hat{\theta}_n-\theta_0)\xrightarrow{d}\mathcal{N}(0,V)\) with \(V=(G^\top WG)^{-1}G^\top WS\,WG(G^\top WG)^{-1}\). The optimal weight \(W=S(\theta_0)^{-1}\) reduces this to the efficient variance \(V_{\mathrm{eff}}=(G^\top S^{-1}G)^{-1}\). 0◻