Robust estimations from distribution structures: II. Central Moments

Tuobang Li a,b,c,1 a Technion-Israel Institute of Technology, Haifa 32000, Israel; b Guangdong Technion-Israel Institute of Technology, Shantou 515063, China; c University of California, Berkeley, CA 94720

This manuscript was compiled on February 17, 2024 In descriptive statistics, U -statistics arise naturally in producing minimum-variance unbiased estimators. In 1984, Serfling considered the distribution formed by evaluating the kernel of the U -statistics and proposed generalized L-statistics which includes Hodges-Lehamnn estimator and Bickel-Lehmann spread as special cases. However, the structures of the kernel distributions remain unclear. In 1954, Hodges and Lehmann demonstrated that if X and Y are independently sampled from the same unimodal distribution, X -Y will exhibit symmetrical unimodality with its peak centered at zero. Building upon this foundational work, the current study delves into the structure of the kernel distribution. It is shown that the kth central moment kernel distributions (k > 2) derived from a unimodal distribution exhibit location invariance and is also nearly unimodal with the mode and median close to zero. This article provides an approach to study the general structure of kernel distributions.

moments | invariant | unimodal | U -statistics | generalized L-statistics
T he most popular robust scale estimator currently, the median absolute deviation, was popularized by [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF] [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF], who credits the idea to Gauss in 1816 [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF]. It can be seen as evaluating the median of a pseudo-sample formed by the absolute deviations of all values related to the sample median. The pseudo-sample size is n. Indeed, most scale estimators can be transformed in such ways. For example, range or interquartile range can be seen as evaluating the mean of a pseudo-sample with two values, and they belong to the class of scale estimators called quantile differences. In 1976, in their landmark series Descriptive Statistics for Nonparametric Models, Bickel and Lehmann (3) generalized a class of estimators as measures of the dispersion of a symmetric distribution around its center of symmetry. Median absolute deviation, sample variance, and average absolute deviation are all belonging to this class. In 1979, the same series, they [START_REF] Bickel | Descriptive statistics for nonparametric models iv[END_REF] proposed a class of estimators referred to as measures of spread, which consider the pairwise differences of a random variable, irrespective of its symmetry, throughout its distribution, rather than focusing on dispersion relative to a fixed point. In the final section (4), they explored a version of the trimmed standard deviation based on n 2 pairwise differences, which is modified here for comparison,

n 1 2 -ϵ -1 2   n(1-ϵ) i= n 2 [Xi -Xn-i+1] 2   1 2
, [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF] and

n 2 (1 -ϵ0 -γϵ0) -1 2   ( n 2 )(1-ϵ0) i=( n 2 )γϵ0 (Xi 1 -Xi 2 ) i 2   1 2
, [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF] where (Xi

1 -Xi 2 ) 1 ≤ . . . ≤ (Xi 1 -Xi 2 ) ( n 2 ) are the order statistics of the pseudosample, Xi 1 -Xi 2 , i1 < i2, provided that n 2 γϵ0 ∈ N and n 2 (1 -ϵ0) ∈ N.
They showed that, when ϵ0 = 0, the result obtained using [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF] is equal to √ 2 times the sample standard deviation. The paper ended with, "We do not know a fortiori which of the measures is preferable and leave these interesting questions open."

Two examples of the impacts of that series are as follows. [START_REF] Oja | On location, scale, skewness and kurtosis of univariate distributions[END_REF][START_REF] Oja | Descriptive statistics for multivariate distributions[END_REF] [START_REF] Oja | On location, scale, skewness and kurtosis of univariate distributions[END_REF]6) provided a more comprehensive and generalized examination of these concepts, and integrated the measures of location, dispersion, and spread as proposed by Bickel and Lehmann [START_REF] Bickel | Descriptive statistics for nonparametric models[END_REF][START_REF] Bickel | Descriptive statistics for nonparametric models iv[END_REF][START_REF] Bickel | Descriptive statistics for nonparametric models ii[END_REF], along with van Zwet's convex transformation order of skewness and kurtosis (1964) [START_REF] Van | Convex transformations: A new approach to skewness and kurtosis in Selected Works of Willem van Zwet[END_REF] for univariate and multivariate distributions, resulting a greater degree of generality and a broader perspective on these statistical constructs. Rousseeuw and Croux proposed a popular efficient scale estimator based on separate medians of pairwise differences taken over i1 and i2 (9) in 1993. However the importance of tackling the symmetry assumption has been greatly underestimated, as will be discussed later.

Here, their open question is addressed in two different aspects [START_REF] Bickel | Descriptive statistics for nonparametric models iv[END_REF]. First, since the estimation of scale can be transformed into the location estimation of a pseudo-sample, according to the principle of the central limit theorem, the variances of such scale estimators should be linearly dependent on the standard deviation of the pseudo-sample and inversely depen-

Significance Statement

In nonparametric statistics, the focus is on the relative differences of robust estimators, which is considered more crucial than their precise values. This principle implies that if the underlying distribution's parameters shift, then all corresponding nonparametric estimates are expected to uniformly and asymptotically adjust in a consistent direction, provided they target the same characteristic of the distribution. This article discusses the validity of this fundamental principle of nonparametrics in various scenarios. It is found that for the kth central moment, kernel distributions generally follow this principle. dent on the square root of the pseudo-sample size. Then, [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF] based on n2 pairwise differences is obviously better than [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF] since the ratio of its pseudo-sample size over that of [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF] is n. So if just considering the size, the variance of [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF] is 1 √ n of the variance of [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF]. Another factor that needs to be considered is the standard deviation of the pseudo-sample. However, the standard deviation of the pseudo-sample is generally independent of different pseudo-sample sizes. So, no matter how different the standard deviation of the pseudo-samples of [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF] and [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF] is, as the sample size increases, the variance of [START_REF] Hampel | The influence curve and its role in robust estimation[END_REF] will always dominate that of [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF]. Second, the nomenclature used in this series is introduced as follows:

Nomenclature. Given a robust estimator, θ, which has an adjustable breakdown point, ϵ, that can approach zero asymptotically, the name of θ comprises two parts: the first part denotes the type of estimator, and the second part represents the population parameter θ, such that θ → θ as ϵ → 0. The abbreviation of the estimator combines the initial letters of the first part and the second part. If the estimator is symmetric, the upper asymptotic breakdown point, ϵ, is indicated in the subscript of the abbreviation of the estimator, with the exception of the median. For an asymmetric estimator based on quantile average, the associated γ follows ϵ.

In REDS I [START_REF] Li | Robust estimations from distribution structures: Mean[END_REF], it was shown that the bias of a robust estimator with an adjustable breakdown point is often monotonic with respect to the breakdown point in a semiparametric distribution. Naturally, the estimator's name should reflect the population parameter that it approaches as ϵ → 0. If multiplying all pseudo-samples by a factor of 1 √ 2 , then [START_REF] Gauss | Bestimmung der genauigkeit der beobachtungen[END_REF] is the trimmed standard deviation adhering to this nomenclature, since ψ2 2 is the kernel function of the unbiased estimation of the second central moment by using U -statistic [START_REF] Heffernan | Unbiased estimation of central moments by using u-statistics[END_REF]. This definition should be preferable, not only because it is the square root of a trimmed U -statistic, which is closely related to the minimum-variance unbiased estimator (MVUE), but also because the second γ-orderliness of the second central moment kernel distribution is ensured by the next exciting theorem.

(x1, x2) = 1 2 (x1 -x2)
Theorem .1. The second central moment kernel distribution generated from any unimodal distribution is second γ-ordered, provided that γ ≥ 0.

Proof. In 1954, Hodges and Lehmann established that if X and Y are independently drawn from the same unimodal distribution, X -Y will be a symmetric unimodal distribution peaking at zero [START_REF] Hodges | Matching in paired comparisons[END_REF]. Given the constraint in the pairwise differences that Xi 1 < Xi 2 , i1 < i2, it directly follows from Theorem 1 in (12) that the pairwise difference distribution (Ξ ∆ ) generated from any unimodal distribution is always monotonic increasing with a mode at zero. Since X -X ′ is a negative variable that is monotonically increasing, applying the squaring transformation, the relationship between the original variable X -X ′ and its squared counterpart (X -X ′ ) 2 can be represented as follows: 2 . In other words, as the negative values of X -X ′ become larger in magnitude (more negative), their squared values (X -X ′ ) 2 become larger as well, but in a monotonically decreasing manner with a mode at zero. Further multiplication by 1 2 also does not change the monotonicity and mode, since the mode is zero. Therefore, the transformed pdf becomes monotonically decreasing with a mode at zero. In REDS I [START_REF] Li | Robust estimations from distribution structures: Mean[END_REF], it was proven that a right-skewed distribution with a monotonic decreasing pdf is always second γ-ordered, which gives the desired result.

X -X ′ < Y -Y ′ =⇒ (X -X ′ ) 2 > (Y -Y ′ )
In REDS I [START_REF] Li | Robust estimations from distribution structures: Mean[END_REF], it was shown that any symmetric distribution is νth U -ordered, suggesting that νth U -orderliness does not require unimodality, e.g., a symmetric bimodal distribution is also νth U -ordered. In the SI Text of REDS I [START_REF] Li | Robust estimations from distribution structures: Mean[END_REF], an analysis of the Weibull distribution showed that unimodality does not assure orderliness. Theorem .1 uncovers a profound relationship between unimodality, monotonicity, and second γ-orderliness, which is sufficient for γ-trimming inequality and γ-orderliness.

On the other hand, while robust estimation of scale has been intensively studied with established methods [START_REF] Bickel | Descriptive statistics for nonparametric models[END_REF][START_REF] Bickel | Descriptive statistics for nonparametric models iv[END_REF], the development of robust measures of asymmetry and kurtosis lags behind, despite the availability of several approaches [START_REF] Al Bowley | Elements of statistics[END_REF][START_REF] Wr Van Zwet | Convex Transformations of Random Variables: Nebst Stellingen[END_REF][START_REF] Ra Groeneveld | Measuring skewness and kurtosis[END_REF][START_REF] Saw | Moments of sample moments of censored samples from a normal population[END_REF][START_REF] Ea Elamir | Trimmed l-moments[END_REF]. The purpose of this paper is to demonstrate that, in light of previous works, the estimation of all central moments can be transformed into a location estimation problem by using U -statistics and the central moment kernel distributions possess desirable properties.

Robust Estimations of the Central Moments

In 1928, Fisher constructed k-statistics as unbiased estimators of cumulants [START_REF] Fisher | Moments and product moments of sampling distributions[END_REF]. Halmos (1946) proved that a functional θ admits an unbiased estimator if and only if it is a regular statistical functional of degree k and showed a relation of symmetry, unbiasness and minimum variance [START_REF] Pr Halmos | The theory of unbiased estimation[END_REF]. Hoeffding, in 1948, generalized U -statistics [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF] which enable the derivation of a minimum-variance unbiased estimator from each unbiased estimator of an estimable parameter. In 1984, Serfling pointed out the speciality of Hodges-Lehmann estimator, which is neither a simple L-statistic nor a U -statistic, and considered the generalized L-statistics and trimmed U -statistics [START_REF] Serfling | Generalized l-, m-, and r-statistics[END_REF]. Given a kernel function h k which is a symmetric function of k variables, the LU -statistic is defined as:

LU h k ,k,k,ϵ,γ,n := LL k,ϵ 0 ,γ,n sort (h k (XN 1 , . . . , XN k )) ( n k ) N =1
,

where ϵ = 1 -(1 -ϵ0) 1 k (proven in Subsection ?? in REDS III (22)), XN 1 , . . . , XN k are the n choose k elements from the sample, LL k,ϵ 0 ,γ,n (Y ) denotes the LL-statistic with the sorted sequence sort (h k (XN 1 , . . . , XN k )) ( n k ) N =1
serving as an input. In the context of Serfling's work, the term 'trimmed U -statistic' is used when LL k,ϵ 0 ,γ,n is TMϵ 0 ,γ,n [START_REF] Serfling | Generalized l-, m-, and r-statistics[END_REF].

In 1997, Heffernan [START_REF] Heffernan | Unbiased estimation of central moments by using u-statistics[END_REF] obtained an unbiased estimator of the kth central moment by using U -statistics and demonstrated that it is the minimum variance unbiased estimator for distributions with the finite first k moments. The weighted [START_REF] Heffernan | Unbiased estimation of central moments by using u-statistics[END_REF]. Despite the complexity, the following theorem offers an approach to infer the general structure of such kernel distributions.

H-L kth central moment (2 ≤ k ≤ n) is thus defined as, WHLkm k,ϵ,γ,n := LU h k =ψ k ,k,k,ϵ,γ,n , where WHLM k,ϵ 0 ,γ,n is used as the LL k,ϵ 0 ,γ,n in LU , ψ k (x1, . . . , x k ) = k-2 j=0 (-1) j 1 k-j x k-j i 1 xi 2 . . . xi j+1 + (-1) k-1 (k -1) x1 . . . x k , the second summation is over i1, . . . , ij+1 = 1 to k with i1 ̸ = i2 ̸ = . . . ̸ = ij+1 and i2 < i3 < . . . < ij+1

Theorem .2. Define a set T comprising all pairs

(ψ k (v), fX,...,X (v)) such that ψ k (v) = ψ k (Q(p1), . . . , Q(p k )) with Q(p1) < . . . < Q(p k ) and fX,...,X (v) = k!f (Q(p1)) . . . f (Q(p k )) is the probability density of the k- tuple, v = (Q(p1), . . . , Q(p k )) (a

formula drawn after a modification of the Jacobian density theorem). T ∆ is a subset of T , consisting all those pairs for which the corresponding

k-tuples satisfy that Q(p1) -Q(p k ) = ∆. The com- ponent quasi-distribution, denoted by ξ ∆ , has a quasi-pdf f ξ ∆ ( ∆) = (ψ k (v),f X,...,X (v))∈T ∆ ∆=ψ k (v) fX,...,X (v), i.e.,
Q(0) -Q(1) to 0. Each component quasi-distribution has a support of - k 3+(-1) k 2 -1 (-∆) k , 1 k (-∆) k .
Proof. The support of ξ ∆ is the extrema of the function

ψ k (Q(p1), ••• , Q(p k )) subjected to the constraints, Q(p1) < ••• < Q(p k ) and ∆ = Q(p1) -Q(p k ).
Using the Lagrange multiplier, the only critical point can be determined at

Q(p1) = ••• = Q(p k ) = 0, where ψ k = 0.
Other candidates are within the boundaries, i.e.,

ψ k (x1 = Q(p1), x2 = Q(p k ), ••• , x k = Q(p k )), •••, ψ k (x1 = Q(p1), ••• , xi = Q(p1), xi+1 = Q(p k ), ••• , x k = Q(p k )), •••, ψ k (x1 = Q(p1), ••• , x k-1 = Q(p1), x k = Q(p k )). ψ k (x1 = Q(p1), ••• , xi = Q(p1), xi+1 = Q(p k ), ••• , x k = Q(p k ))
can be divided into k groups. The gth group has the common factor (-1) g+1 1 k-g+1 , if 1 ≤ g ≤ k -1 and the final kth group is the term (-1)

k-1 (k -1) Q(p1) i Q(p k ) k-i . If k+1-i 2 ≤ j ≤ k-1 2 and j + 1 ≤ g ≤ k -j, the gth group has i i-1 g-j-1 k-i j terms having the form (-1) g+1 1 k-g+1 Q(p1) k-j Q(p k ) j . If k+1-i 2 ≤ j ≤ k-1 2 and k -j + 1 ≤ g ≤ i + j, the gth group has i i-1 g-j-1 k-i j + (k -i) k-i-1 j-k+g-1 i k-j terms having the form (-1) g+1 1 k-g+1 Q(p1) k-j Q(p k ) j . If 0 ≤ j < k+1-i 2 and j+1 ≤ g ≤ i+j, the gth group has i i-1 g-j-1 k-i j terms having the form (-1) g+1 1 k-g+1 Q(p1) k-j Q(p k ) j . If k 2 ≤ j ≤ k and k -j + 1 ≤ g ≤ j, the gth group has (k -i) k-i-1 j-k+g-1 i k-j terms having the form (-1) g+1 1 k-g+1 Q(p1) k-j Q(p k ) j . If k 2 ≤ j ≤ k and j + 1 ≤ g ≤ j + i < k, the gth group has i i-1 g-j-1 k-i j +(k -i) k-i-1 j-k+g-1 i k-j terms having the form (-1) g+1 1 k-g+1 Q(p1) k-j Q(p k ) j . So, if i + j = k, k 2 ≤ j ≤ k, 0 ≤ i ≤ k 2 , the summed coefficient of Q(p1) i Q(p k ) k-i is (-1) k-1 (k -1) + k-1 g=i+1 (-1) g+1 1 k-g+1 (k -i) k-i-1 g-i-1 + k-1 g=k-i+1 (-1) g+1 1 k-g+1 i i-1 g-k+i-1 = (-1) k-1 (k -1) + (-1) k+1 + (k -i) (-1) k + (-1) k (i -1) = (-1) k+1 .
The summation identities are

k-1 g=i+1 (-1) g+1 1 k-g+1 (k -i) k-i-1 g-i-1 = (k -i) 1 0 k-1 g=i+1 (-1) g+1 k-i-1 g-i-1 t k-g dt = (k -i) 1 0 (-1) i (t -1) k-i-1 -(-1) k+1 dt = (k -i) (-1) k i-k + (-1) k = (-1) k+1 + (k -i) (-1) k and k-1 g=k-i+1 (-1) g+1 1 k-g+1 i i-1 g-k+i-1 = 1 0 k-1 g=k-i+1 (-1) g+1 i i-1 g-k+i-1 t k-g dt = 1 0 i (-1) k-i (t -1) i-1 -i (-1) k+1 dt = (-1) k (i -1). If 0 ≤ j < k+1-i 2 and i = k, ψ k = 0. If k+1-i 2 ≤ j ≤ k-1 2 and k+1 2 ≤ i ≤ k -1, the summed coefficient of Q(p1) i Q(p k ) k-i is (-1) k-1 (k -1) + k-1 g=k-i+1 (-1) g+1 1 k-g+1 i i-1 g-k+i-1 + k-1 g=i+1 (-1) g+1 1 k-g+1 (k -i) k-i-1 g-i-1 , the same as above. If i + j < k, since i k-j
= 0, the related terms can be ignored, so, using the binomial theorem and beta function, the summed coefficient of

Q(p1) k-j Q(p k ) j is i+j g=j+1 (-1) g+1 1 k-g+1 i i-1 g-j-1 k-i j = i k-i j 1 0 i+j g=j+1 (-1) g+1 i-1 g-j-1 t k-g dt = k-i j i 1 0 (-1) j t k-j-1 t t-1 1-i dt = k-i j i (-1) j+i+1 Γ(i)Γ(k-j-i+1) Γ(k-j+1) = (-1) j+i+1 i!(k-j-i)!(k-i)! (k-j)!j!(k-j-i)! = (-1) j+i+1 i!(k-i)! k! k! (k-j)!j! = k i -1 (-1) 1+i k j (-1) j . According to the binomial theorem, the coefficient of Q(p1) i Q(p k ) k-i in k i -1 (-1) 1+i (Q(p1) -Q(p k )) k is k i -1 (-1) 1+i k i (-1) k-i = (-1) k+1 , same as the above summed coefficient of Q(p1) i Q(p k ) k-i , if i + j = k. If i + j < k, the coefficient of Q(p1) k-j Q(p k ) j is k i -1 (-1) 1+i k j (-1) j , same as the corresponding summed coefficient of Q(p1) k-j Q(p k ) j .
Therefore,

ψ k (x1 = Q(p1), . . . , xi = Q(p1), xi+1 = Q(p k ), . . . , x k = Q(p k )) = k i -1 (-1) 1+i (Q(p1) -Q(p k )) k ,
the maximum and minimum of ψ k follow directly from the properties of the binomial coefficient.

The component quasi-distribution, ξ ∆ , is closely related to Ξ ∆ , which is the pairwise difference distribution, since

1 k (-∆) k ∆=-( k 3+(-1) k 2 ) -1 (-∆) k f ξ ∆ ( ∆) = f Ξ ∆ (∆).
Recall that Theorem .1 established that f Ξ ∆ (∆) is monotonic increasing with a mode at zero if the original distribution is unimodal, f Ξ -∆ (-∆) is thus monotonic decreasing with a mode at zero. In general, if assuming the shape of ξ ∆ is uniform, Ξ k is monotonic left and right around zero. The median of Ξ k also exhibits a strong tendency to be close to zero, as it can be cast as a weighted mean of the medians of ξ ∆ . When -∆ is small, all values of ξ ∆ are close to zero, resulting in the median of ξ ∆ being close to zero as well. When -∆ is large, the median of ξ ∆ depends on its skewness, but the corresponding weight is much smaller, so even if ξ ∆ is highly skewed, the median of Ξ k will only be slightly shifted from zero. Denote the median of Ξ k as mkm, for the five parametric distributions here, |mkm|s are all ≤ 0.1σ for Ξ3 and Ξ4, where σ is the standard deviation of Ξ k (SI Dataset S1). Assuming mkm = 0, for the even ordinal central moment kernel distribution, the average probability density on the left side of zero is greater than that on the right side, since

1 2 ( k 2 ) -1 (Q(0)-Q(1)) k > 1 2 1 k (Q(0)-Q(1)
) k . This means that, on average, the inequality f (Q(ϵ)) ≥ f (Q(1 -ϵ)) holds. For the odd ordinal distribution, the discussion is more challenging since it is generally symmetric. Just consider Ξ3, let x1 = Q(pi) and x3 = Q(pj), changing the value of x2 from Q(pi) to Q(pj) will monotonically change the value of ψ3(x1, x2, x3),

since ∂ψ 3 (x 1 ,x 2 ,x 3 ) ∂x 2 = - x 2 1 2 -x1x2 + 2x1x3 + x 2 2 -x2x3 - x 2 3 2 , -3 4 (x1 -x3) 2 ≤ ∂ψ 3 (x 1 ,x 2 ,x 3 ) ∂x 2 ≤ -1 2 (x1 -x3) 2 ≤ 0.
If the original distribution is right-skewed, ξ ∆ will be left-skewed, so, for Ξ3, the average probability density of the right side of zero will be greater than that of the left side, which means, on average, the inequality f (Q(ϵ)) ≤ f (Q(1 -ϵ)) holds. In all, the monotonic decreasing of the negative pairwise difference distribution guides the general shape of the kth central moment kernel distribution, k > 2, forcing it to be unimodal-like with the mode and median close to zero, then, the inequality

f (Q(ϵ)) ≤ f (Q(1 -ϵ)) or f (Q(ϵ)) ≥ f (Q(1 -ϵ)) holds in general.
If a distribution is νth γ-ordered and all of its central moment kernel distributions are also νth γ-ordered, it is called completely νth γ-ordered.

Another crucial property of the central moment kernel distribution, location invariant, is introduced in the next theorem.

Theorem .3. ψ k (x1 = λx1 + µ, ••• , x k = λx k + µ) = λ k ψ k (x1, ••• , x k ).
Proof. Recall that for the kth central moment, the kernel is

ψ k (x1, . . . , x k ) = k-2 j=0 (-1) j 1 k-j x k-j i 1 xi 2 . . . xi j+1 + (-1) k-1 (k -1) x1 . . . x k ,
where the second summation is over i1, . . . , ij+1 = 1 to k with i1 ̸ = i2 ̸ = . . . ̸ = ij+1 and i2 < i3 < . . . < ij+1 [START_REF] Heffernan | Unbiased estimation of central moments by using u-statistics[END_REF].

ψ k consists of two parts. The first part,

k-2 j=0 (-1) j 1 k-j
x k-j i 1 xi 2 . . . xi j+1 , involves a double summation over certain terms.

The second part, (-1) k-1 (k -1) x1 . . . x k , carries an alternating sign (-1) k-1 and involves multiplication of the constant k -1 with the product of all the x variables, x1x2 . . . x k . Consider each multiplication cluster (-1) j 1 k-j

x k-j i 1 xi 2 . . . xi j+1 for j ranging from 0 to k -2 in the first part. Let each cluster form a single group. The first part can be divided into k -1 groups. Combine this with the second part (-1) k-1 (k -1) x1 . . . x k . Together, the terms of ψ k can be divided into a total of k groups. From the 1st to k -1th group, the gth group has k g g 1 terms having the form

(-1) g+1 1 k-g+1 x k-g+1 i 1 xi 2 . . . xi g . The final kth group is the term (-1) k-1 (k -1) x1 ••• x k .
There are two ways to divide ψ k into k groups according to the form of each term. The first choice is, if k ̸ = g, the gth group of ψ k has k-l g-l terms having the form (-1)

g+1 1 k-g+1 x k-g+1 i 1 xi 2 ••• xi l xi l+1 . . . xi g , where xi 1 , xi 2 , ••• , xi l are fixed, xi l+1 , ••• , xi g are selected such that i l+1 , ••• , ig ̸ = i1, i2, ••• , i l and i l+1 ̸ = . . . ̸ = ig. De- fine another function Ψ k xi 1 , xi 2 , ••• , xi l , xi l+1 , ••• , xi g = (λxi 1 + µ) k-g+1 (λxi 2 + µ) ••• (λxi l + µ) λxi l+1 + µ ••• λxi g + µ , the first group of Ψ k is λ k xi 1 ••• xi l xi l+1 ••• xi g , the hth group of Ψ k , h > 1, has k-g+1 k-h-l+2 terms hav- ing the form λ k-h+1 µ h-1 x k-h-l+2 i 1 xi 2 ••• xi l .
Transforming ψ k by Ψ k , then combing all terms with

λ k-h+1 µ h-1 x k-h-l+2 i 1 xi 2 ••• xi l , k -h -l + 2 > 1, the summed coefficient is S1 l = h+l-1 g=l (-1) g+1 1 k-g+1 k-g+1 k-h-l+2 k-l g-l = h+l-1 g=l (-1) g+1 (k-l)! (h+l-g-1)!(k-h-l+2)!(g-l)!
= 0, since the summation is starting from l, ending at h + l -1, the first term includes the factor g -l = 0, the final term includes the factor h + l -g -1 = 0, the terms in the middle are also zero due to the factorial property.

Another possible choice is the gth group of

ψ k has (k -h) h-1 g-k+h-1 terms having the form (-1) g+1 1 k-g+1 xi 1 xi 2 ••• x k-g+1 i j ••• xi k-h+1 xi k-h+2 ••• xi g , provided that k ̸ = g, 2 ≤ j ≤ k -h + 1, where xi 1 , . . . , xi k-h+1 are fixed, x k-g+1 i j and xi k-h+2 , ••• , xi g are selected such that i k-h+2 , ••• , ig ̸ = i1, i2, ••• , i k-h+1 and i k-h+2 ̸ = . . . ̸ = ig. Transforming these terms by Ψ k xi 1 , xi 2 , . . . , xi j , . . . , xi k-h+1 , xi k-h+2 , . . . , xi g = (λxi 1 + µ) (λxi 2 + µ) ••• λxi j + µ k-g+1 ••• λxi k-h+1 + µ λxi k-h+2 + µ
then there are k -g + 1 terms having the form

λ k-h+1 µ h-1 xi 1 xi 2 . . . xi k-h+1 . Transforming the final kth group of ψ k by Ψ k (x1, . . . , x k ) = (λx1 + µ) ••• (λx k + µ), then, there is one term having the form (-1) k-1 (k -1) λ k-h+1 µ h-1 x1x2 . . . x k-h+1 . An- other possible combination is that the gth group of ψ k contains (g -k + h -1) h-1 g-k+h-1
terms having the form

(-1) g+1 1 k-g+1 xi 1 xi 2 ••• xi k-h+1 xi k-h+2 ••• x k-g+1 i j ••• xi g . Transforming these terms by Ψ k xi 1 , xi 2 , . . . , xi k-h+1 , xi k-h+2 , . . . , xi j , . . . , xi g = (λxi 1 + µ) (λxi 2 + µ) ••• λxi k-h+1 + µ λxi k-h+2 + µ ••• λxi j + µ k-g+1
then there is only one term having the form

λ k-h+1 µ h-1 xi 1 xi 2 . . . xi k-h+1 .
The above summation S1 l should also be included, i.e., x k-h-l+2

i 1 = xi 1 , k = h + l -1. So, combing all terms with λ k-h+1 µ h-1 xi 1 xi 2 . . . xi k-h+1 , accord- ing to the binomial theorem, the summed coefficient is S2 l = k-1 g=k-h+1 (-1) g+1 h-1 g-k+h-1 k -h + 1 + g-k+h-1 k-g+1 + (-1) k-1 (k -1) = (k -h + 1) k-1 g=k-h+1 (-1) g+1 h-1 g-k+h-1 + k-1 g=k-h+1 (-1) g+1 h-1 g-k+h-1 g-k+h-1 k-g+1 + (-1) k-1 (k -1) = (-1) k (k -h + 1) + (h -2)(-1) k + (-1) k-1 (k -1) = 0.
The summation identities required are

k-1 g=k-h+1 (-1) g+1 h-1 g-k+h-1 = (-1) k and k-1 g=k-h+1 (-1) g+1 h-1 g-k+h-1 g-k+h-1 k-g+1
= (h -2)(-1) k . These two summation identities are proven in Lemma 4 and 5 in the SI Text.

Thus, no matter in which way, all terms including µ can be canceled out. The proof is complete by noticing that the remaining part is

λ k ψ k (x1, ••• , x k ).
A direct result of Theorem .3 is that, WHLkm after standardization is invariant to location and scale. So, the weighted H-L standardized kth moment is defined to be

WHLskm ϵ=min (ϵ 1 ,ϵ 2 ),k 1 ,k 2 ,γ 1 ,γ 2 ,n := WHLkm k 1 ,ϵ 1 ,γ 1 ,n (WHLvar k 2 ,ϵ 2 ,γ 2 ,n ) k/2 .
To avoid confusion, it should be noted that the robust location estimations of the kernel distributions discussed in this paper differ from the approach taken by Joly and Lugosi (2016) [START_REF] Joly | Robust estimation of u-statistics[END_REF], which is computing the median of all U -statistics from different disjoint blocks. Compared to bootstrap median U -statistics, this approach can produce two additional kinds of finite sample bias, one arises from the limited numbers of blocks, another is due to the size of the U -statistics (consider the mean of all U -statistics from different disjoint blocks, it is definitely not identical to the original U -statistic, except when the kernel is the Hodges-Lehmann kernel). Laforgue, Clemencon, and Bertail (2019)'s median of randomized Ustatistics ( 24) is more sophisticated and can overcome the limitation of the number of blocks, but the second kind of bias remains unsolved.

Congruent Distribution

In the realm of nonparametric statistics, the relative differences, or orders, of robust estimators are of primary importance. A key implication of this principle is that when there is a shift in the parameters of the underlying distribution, all nonparametric estimates should asymptotically change in the same direction, if they are estimating the same attribute of the distribution. If, on the other hand, the mean suggests an increase in the location of the distribution while the median indicates a decrease, a contradiction arises. It is worth noting that such contradiction is not possible for any LL-statistics in a location-scale distribution, as explained in Theorem 2 and 18 in REDS I. However, it is possible to construct counterexamples to the aforementioned implication in a shape-scale distribution. In the case of the Weibull distribution, its quantile function is QW ei (p) = λ(-ln(1 -p)) 1/α , where 0 ≤ p ≤ 1, α > 0, λ > 0, λ is a scale parameter, α is a shape parameter, ln is the natural logarithm function. Then, m = λ α ln(2), µ = λΓ 1 + 1 α , where Γ is the gamma function. When α = 1, m = λ ln(2) ≈ 0.693λ, µ = λ, when α = 1 2 , m = λ ln 2 (2) ≈ 0.480λ, µ = 2λ, the mean increases as α changes from 1 to 1 2 , but the median decreases. In the last section, the fundamental role of quantile average was demonstrated by using the method of classifying distributions through the signs of derivatives. To avoid such scenarios, this method can also be used. Let the quantile average function of a parametric distribution be denoted as QA (ϵ, γ, α1,

••• , αi, ••• , α k )
, where αi represent the parameters of the distribution, then, a distribution is γ-congruent if and only if the sign of ∂QA ∂α i remains the same for all 0 ≤ ϵ ≤ 1 1+γ . If ∂QA ∂α i is equal to zero or undefined, it can be considered both positive and negative, and thus does not impact the analysis. A distribution is completely γ-congruent if and only if it is γ-congruent and all its central moment kernel distributions are also γ-congruent. Setting γ = 1 constitutes the definitions of congruence and complete congruence. Replacing the QA with QHLM gives the definition of γ-U -congruence. Chebyshev's inequality implies that, for any probability distributions with finite second moments, as the parameters change, even if some LL-statistics change in a direction different from that of the population mean, the magnitude of the changes in the LL-statistics remains bounded compared to the changes in the population mean. Furthermore, distributions with infinite moments can be γ-congruent, since the definition is based on the quantile average, not the population mean.

The following theorems show the conditions that a distribution is congruent or γ-congruent.

Theorem .4. A symmetric distribution is always congruent and U -congruent.

Proof. As shown in Theorem 2 and Theorem 18 in REDS I, for any symmetric distribution, all symmetric quantile averages and all SQHLMs conincide. The conclusion follows immediately.

Theorem .5. A positive definite location-scale distribution is always γ-congruent.

Proof. As shown in Theorem 2, for a location-scale distribution, any quantile average can be expressed as λQA 0 (ϵ, γ) + µ. Therefore, the derivatives with respect to the parameters λ or µ are always positive. By application of the definition, the desired outcome is obtained.

For the Pareto distribution, ∂Q ∂α = xm(1-p) -1/α ln(1-p) α 2

. Since ln(1 -p) < 0 for all 0 < p < 1, (1 -p) -1/α > 0 for all 0 < p < 1 and α > 0, so ∂Q ∂α < 0, and therefore ∂QA ∂α < 0, the Pareto distribution is γcongruent. It is also γ-U -congruent, since γmHLM can also express as a function of Q(p). For the lognormal dis-

tribution, ∂QA ∂σ = 1 2 √ 2erfc -1 (2γϵ) -e √ 2µ-2σerfc -1 (2γϵ) √ 2 + - √ 2 erfc -1 (2(1 -ϵ))e √ 2µ-2σerfc -1 (2(1-ϵ)) √ 2
. Since the inverse complementary error function is positive when the input is smaller than 1, and negative when the input is larger than 1, and symmetry around 1

, if 0 ≤ γ ≤ 1, erfc -1 (2γϵ) ≥ -erfc -1 (2 -2ϵ), e µ- √ 2σerfc -1 (2-2ϵ) > e µ- √ 2σerfc -1 (2γϵ) . Therefore, if 0 ≤ γ ≤ 1, ∂QA ∂σ
> 0, the lognormal distribution is γ-congruent. Theorem .4 implies that the generalized Gaussian distribution is congruent and U -congruent. For the Weibull distribution, when α changes from 1 to 1 2 , the average probability density on the left side of the median increases, since

1 2 λ ln(2) < 1 2
λ ln 2 (2) , but the mean increases, indicating that the distribution is more heavy-tailed, the probability density of large values will also increase. So, the reason for non-congruence of the Weibull distribution lies in the simultaneous increase of probability densities on two opposite sides as the shape parameter changes: one approaching the bound zero and the other approaching infinity. Note that the gamma distribution does not have this issue, Numerical results indicate that it is likely to be congruent.

The next theorem shows an interesting relation between congruence and the central moment kernel distribution.

Theorem .6. The second central moment kernal distribution derived from a continuous location-scale unimodal distribution is always γ-congruent.

Proof. Theorem .3 shows that the central moment kernel distribution generated from a location-scale distribution is also a location-scale distribution. Theorem .1 shows that it is positively definite. Implementing Theorem 12 in REDS 1 yields the desired result.

Although some parametric distributions are not congruent, as shown in REDS 1. In REDS 1, Theorem 12 establishes that γ-congruence always holds for a positive definite location-scale family distribution and thus for the second central moment kernel distribution generated from a location-scale unimodal distribution as shown in Theorem .6. Theorem .2 demonstrates that all central moment kernel distributions are unimodal-like with mode and median close to zero, as long as they are generated from unimodal distributions. Assuming finite moments and constant Q(0)-Q(1), increasing the mean of a distribution will result in a generally more heavy-tailed distribution, i.e., the probability density of the values close to Q(1) increases, since the total probability density is 1. In the case of the kth central moment kernel distribution, k > 2, while the total probability density on either side of zero remains generally constant as the median is generally close to zero and much less impacted by increasing the mean, the probability density of the values close to zero decreases as the mean increases. This transformation will increase nearly all symmetric weighted averages, in the general sense. Therefore, except for the median, which is assumed to be zero, nearly all symmetric weighted averages for all central moment kernel distributions derived from unimodal distributions should change in the same direction when the parameters change.

Discussion

Moments, including raw moments, central moments, and standardized moments, are the most common parameters that describe probability distributions. Central moments are preferred over raw moments because they are invariant to translation. In 1947, Hsu and Robbins proved that the arithmetic mean converges completely to the population mean provided the second moment is finite [START_REF] Hsu | Complete convergence and the law of large numbers[END_REF]. The strong law of large numbers (proven by Kolmogorov in 1933) [START_REF] Kolmogorov | Sulla determinazione empirica di una lgge di distribuzione[END_REF] implies that the kth sample central moment is asymptotically unbiased. Recently, fascinating statistical phenomena regarding Taylor's law for distributions with infinite moments have been discovered by [START_REF] Drton | Wald tests of singular hypotheses[END_REF] [START_REF] Drton | Wald tests of singular hypotheses[END_REF], Pillai and Meng (2016) [START_REF] Ns Pillai | An unexpected encounter with cauchy and lévy[END_REF], Cohen, Davis, and Samorodnitsky (2020) [START_REF] Je Cohen | Heavy-tailed distributions, correlations, kurtosis and taylor's law of fluctuation scaling[END_REF], and Brown, Cohen, Tang, and Yam (2021) [START_REF] Brown | Taylor's law of fluctuation scaling for semivariances and higher moments of heavy-tailed data[END_REF]. Lindquist and Rachev (2021) raised a critical question in their inspiring comment to Brown et al's paper [START_REF] Brown | Taylor's law of fluctuation scaling for semivariances and higher moments of heavy-tailed data[END_REF]: "What are the proper measures for the location, spread, asymmetry, and dependence (association) for random samples with infinite mean?" [START_REF] Lindquist | Taylor's law and heavy-tailed distributions[END_REF]. From a different perspective, this question closely aligns with the essence of Bickel and Lehmann's open question in 1979 (4). They suggested using median, interquartile range, and medcouple [START_REF] Brys | A robust measure of skewness[END_REF] as the robust versions of the first three moments. While answering this question is not the focus of this paper, it is almost certain that the estimators proposed in this paper will have a place. Since the estimation of central moments can be transformed into the location estimation of a pseudosample, according to the general principle of central limit theorem, the optimal estimator should always has a combinatorial pseudosample size, which explains, in another aspect, why the theory of U -statistics allows a minimum-variance unbiased estimator to be derived from each unbiased estimator of an estimable parameter. Similar to the robust version of Lmoment [START_REF] Hosking | Analysis and estimation of distributions using linear combinations of order statistics[END_REF] being trimmed L-moment [START_REF] Ea Elamir | Trimmed l-moments[END_REF], central moments now also have their robust nonparametric version, weighted Hodges-Lehmann central moments, based on the complete U -congruence of the underlying distribution.

Software Availability. The codes used to compute the weighted H-L kth central moment have been deposited in GitHub.

To whom correspondence should be addressed. E-mail: tl@biomathematics.org PNAS | February 17,

| vol. XXX | no. XX | 1-6

| Li

| Li


Robust estimations from distribution structures: III. Invariant Moments

Tuobang Li a,b,c,1 a Technion-Israel Institute of Technology, Haifa 32000, Israel; b Guangdong Technion-Israel Institute of Technology, Shantou 515063, China; c University of California, Berkeley, CA 94720

This manuscript was compiled on [START_REF]Reproducibility and Replicability in Science[END_REF] Descriptive statistics for parametric models are currently highly sensative to departures, gross errors, and/or random errors. Here, leveraging the structures of parametric distributions and their central moment kernel distributions, a class of estimators, consistent simultanously for both a semiparametric distribution and a distinct parametric distribution, is proposed. These efficient estimators are robust to both gross errors and departures from parametric assumptions, making them ideal for estimating the mean and central moments of common unimodal distributions. This article also illuminates the understanding of the common nature of probability distributions and the measures of them.

moments | invariant | unimodal | adaptive estimation | U -statistics
T he potential biases of robust location estimators in esti- mating the population mean have been noticed for more than two centuries [START_REF] Gauss | Theoria combinationis observationum erroribus minimis obnoxiae[END_REF], with numerous significant attempts made to address them. In calculating a robust estimator, the procedure of identifying and downweighting extreme values inherently necessitates the formulation of distributional assumptions. Previously, it was demonstrated that, due to the presence of infinite-dimensional nuisance shape parameters, the semiparametric approach struggles to consistently address distributions with shapes more intricate than symmetry. Newcomb (1886) provided the first modern approach to robust parametric estimation by developing a class of estimators that gives "less weight to the more discordant observations" [START_REF] Newcomb | A generalized theory of the combination of observations so as to obtain the best result[END_REF]. In 1964, Huber (3) used the minimax procedure to obtain M -estimator for the contaminated normal distribution, which has played a pre-eminent role in the later development of robust statistics. However, as previously demonstrated, under growing asymmetric departures from normality, the bias of the Huber M -estimator (HM) increases rapidly. This is a common issue in parametric robust statistics. For example, He and Fung (1999) constructed (4) a robust M -estimator (HFM) for the two-parameter Weibull distribution, from which the mean and central moments can be calculated. Nonetheless, it is inadequate for other parametric distributions, e.g., the gamma, Perato, lognormal, and the generalized Gaussian distributions (SI Dataset S1). Another interesting approach is based on L-estimators, such as percentile estimators. For examples of percentile estimators for the Weibull distribution, the reader is referred to the works of [START_REF] Menon | Estimation of the shape and scale parameters of the weibull distribution[END_REF] [START_REF] Menon | Estimation of the shape and scale parameters of the weibull distribution[END_REF], [START_REF] Sd Dubey | Some percentile estimators for weibull parameters[END_REF] [START_REF] Sd Dubey | Some percentile estimators for weibull parameters[END_REF], [START_REF] Marks | Estimation of weibull parameters from common percentiles[END_REF] [START_REF] Marks | Estimation of weibull parameters from common percentiles[END_REF], and Boudt, Caliskan, and Croux (2011) [START_REF] Boudt | Robust explicit estimators of weibull parameters[END_REF]. At the outset of the study of percentile estimators, it was known that they arithmetically utilize the invariant structures of parametric distributions [START_REF] Menon | Estimation of the shape and scale parameters of the weibull distribution[END_REF][START_REF] Sd Dubey | Some percentile estimators for weibull parameters[END_REF]). An estimator is classified as an I-statistic if it asymptotically satisfies I (LE1, . . . , LE l ) = (θ1, . . . , θq) for the distribution it is consistent, where LEs are calculated with the use of LU -statistics (defined in REDS II), I is defined using arithmetic operations and constants but may also incorporate transcendental functions and quantile functions, and θs are the population parameters it estimates. In this article, two subclasses of Istatistics are introduced, recombined I-statistics and quantile I-statistics. Based on LU -statistics, I-statistics are naturally robust. Compared to probability density functions (pdfs) and cumulative distribution functions (cdfs), the quantile functions of many parametric distributions are more elegant. Since the expectation of an L-estimator can be expressed as an integral of the quantile function, I-statistics are often analytically obtainable. However, it is observed that even when the sample follows a gamma distribution, which belongs to the same larger family as the Weibull model, the generalized gamma distribution, a misassumption can still lead to substantial biases in Marks percentile estimator (MP) for the Weibull distribution [START_REF] Marks | Estimation of weibull parameters from common percentiles[END_REF] (SI Dataset S1).

The purpose of this paper is to demonstrate that, in light of previous works, by utilizing the invariant structures of unimodal distributions, a suite of robust estimators can be constructed whose biases are typically smaller than the variances (as seen in Table 1 for n = 5184).

A. Invariant Moments. Most popular robust location estimators, such as the symmetric trimmed mean, symmetric Winsorized mean, Hodges-Lehmann estimator, Huber Mestimator, and median of means, are symmetric. As shown in REDS I, a symmetric weighted Hodges-Lehmann mean (SWHLM k,ϵ ) can achieve consistency for the population mean in any symmetric distribution with a finite mean. However, it falls considerably short of consistently handling other parametric distributions that are not symmetric. Shifting from semiparametrics to parametrics, consider a robust estimator with a non-sample-dependent breakdown point (defined in Subsection F) which is consistent simultaneously for both a semiparametric distribution and a parametric distribution that does not belong to that semiparametric distribution, it

Significance Statement

Bias, variance, and contamination are the three main errors in statistics. Consistent robust estimation is unattainable without distributional assumptions. In this article, invariant moments are proposed as a means of achieving near-consistent and robust estimations of moments, even in scenarios where moderate violations of distributional assumptions occur, while the variances are sometimes smaller than those of the sample moments.

T.L. designed research, performed research, analyzed data, and wrote the paper.

The author declares no competing interest.

is named with the prefix 'invariant' followed by the name of the population parameter it is consistent with. Here, the recombined I-statistic is defined as

RI d,h k ,k 1 ,k 2 ,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,n,LU 1 ,LU 2 := lim c→∞ LU1 h k ,k 1 ,k 1 ,ϵ 1 ,γ 1 ,n + c d+1 LU2 h k ,k 2 ,k 2 ,ϵ 2 ,γ 2 ,n + c d -c ,
where d is the key factor for bias correction, LU h k ,k,k,ϵ,γ,n is the LU -statistic, k is the degree of the U -statistic, k is the degree of the LL-statistic, ϵ is the upper asymptotic breakdown point of the LU -statistic. It is assumed in this series that in the subscript of an estimator, if k, k and γ are omitted, First consider the location-scale distribution.

k = 1, k = 1, γ = 1 are assumed, if just one k is indicated, k1 = k2, if just one γ is indicated, γ1 = γ2, if n is omitted,
rm d,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,n,WL 1 ,WL 2 := RI d,h k =x,k 1 =1,k 2 =1,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,n,LU 1 =WL 1 ,LU 2 =WL 2 . Assuming finite means, rm d= µ-WL 1k 1 ,ϵ 1 ,γ 1 WL 1 k 1 ,ϵ 1 ,γ 1 -WL 2 k 2 ,ϵ 2 ,γ 2 ,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,WL
Since

rm d,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,WL 1 ,WL 2 = limc→∞ (WL1k 1 ,ϵ 1 ,γ 1 +c) d+1 (WL2k 2 ,ϵ 2 ,γ 2 +c) d -c = (d + 1) WL1 k 1 ,ϵ 1 ,γ - dWL2 k 2 ,ϵ 2 ,γ = µ. So, d = µ-WL 1k 1 ,ϵ 1 ,γ 1 WL 1k 1 ,ϵ 1 ,γ 1 -WL 2k 2 ,ϵ 2 ,γ 2
. In REDS I, it was established that any WL(k, ϵ, γ) can be expressed as λWL0(k, ϵ, γ) + µ for a location-scale distribution parameterized by a location parameter µ and a scale parameter λ, where WL0(k, ϵ, γ) is a function of Q0(p), the quantile function of a standard distribution without any shifts or scaling, according to the definition of the weighted L-statistic.

The simultaneous cancellation of µ and λ in

(λµ 0 +µ)-(λWL 1 0 (k 1 ,ϵ 1 ,γ 1 )+µ) (λWL 10 (k 1 ,ϵ 1 ,γ 1 )+µ)-(λWL 20 (k 2 ,ϵ 2 ,γ 2 )+µ)
assures that the d in rm is always a constant for a location-scale distribution. The proof of the second assertion follows directly from the coincidence property.

According to Theorem 17 in REDS I, for any symmetric distribution with a finite mean, SWHLM1

k 1 = SWHLM2 k 2 = µ. Then rm d,k 1 ,k 2 ,ϵ 1 ,ϵ 2 ,SWHLM 1 ,SWHLM 2 = limc→∞ (µ+c) d+1
(µ+c) d -c = µ. This completes the demonstration.

For example, the Pareto distribution has a quantile function

QP ar (p) = xm(1 -p) -1
α , where xm is the minimum possible value that a random variable following the Pareto distribution can take, serving a scale parameter, α is a shape parameter. The mean of the Pareto distribution is given by αxm α-1 . As WL(k, ϵ, γ) can be expressed as a function of Q(p), one can set the two WL k,ϵ,γ s in the d value of rm as two arbitrary quantiles QP ar (p1) and QP ar (p2). For the Pareto distribution, dP er,rm = µ P er -Q P ar (p 1 )

Q P ar (p 1 )-Q P ar (p 2 ) = αxm α-1 -xm(1-p 1 ) -1 α xm(1-p 1 ) -1 α -xm(1-p 2 ) -1 α .
xm can be canceled out. Intriguingly, the quantile function of exponential distribution is Qexp(p) = ln 1 1-p λ, λ ≥ 0. µexp = λ. Then, dexp,rm = µexp-Qexp(p 1 ) ,BM,m for distributions with skewness between those of the exponential and symmetric distributions are tiny (SI Dataset S1). rm d≈0.103,ν=3,ϵ= 1 24 ,BM,m exhibits excellent performance for all these common unimodal distributions (SI Dataset S1).

Qexp(p 1 )-Qexp(p 2 ) = λ-ln 1 1-p 1 λ ln 1 1-p 1 λ-ln 1 1-p 2 λ = - ln(1-p 1 )+1 ln(1-p 1 )-ln(1-p 2 ) . Since limα→∞ α α-1 -(1-p 1 ) -1/α (1-p 1 ) -1/α -(1-p 2 ) -1/α = - ln(1-p 1 )+1
d= µ-SWHLM 1k 1 ,ϵ 1 SWHLM 1k 1 ,ϵ 1 -SWHLM 2 k 2 ,ϵ 2 ,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),SWHLM
The recombined mean is a recombined I-statistic. Consider an I-statistic whose LEs are percentiles of a distribution obtained by plugging LU -statistics into a cumulative distribution function, I is defined with arithmetic operations, constants, and quantile functions, such an estimator is classified as a quantile I-statistic. One version of the quantile I-statistic can be defined as QI d,h k ,k,k,ϵ,γ,,n,LU :

= Qn,h k Fn,h k (LU ) -γ 1+γ d + Fn,h k (LU ) Fn,h k (LU ) ≥ γ 1+γ Qn,h k Fn,h k (LU ) -γ 1+γ -Fn,h k (LU ) d Fn,h k (LU ) < γ 1+γ , where LU is LU k,k,ϵ,γ,n , Fn,h k (x) is the empirical cumulative distribution function of the h k kernel distribution, Qn,h k is the quantile function of the h k kernel distribution.
Similarly, the quantile mean can be defined as qm d,k,ϵ,γ,n,WL := QI d,h k =x,k=1,k,ϵ,γ,n,LU =WL . Moreover, in extreme right-skewed heavy-tailed distributions, if the calculated percentile exceeds 1 -ϵ, it will be adjusted to 1 -ϵ. In a left-skewed distribution, if the obtained percentile is smaller than γϵ, it will also be adjusted to γϵ. Without loss of generality, in the following discussion, only the case where Fn (WL k,ϵ,γ,n ) ≥ γ 1+γ is considered. The most popular method for computing the sample quantile function was proposed by Hyndman and Fan in 1996 [START_REF] Rj Hyndman | Sample quantiles in statistical packages[END_REF]. Another widely used method for calculating the sample quantile function in-volves employing linear interpolation of modes corresponding to the order statistics of the uniform distribution on the interval [0, 1], i.e., Qn (p) = X ⌊h⌋ + (h -⌊h⌋) X ⌈h⌉ -X ⌊h⌋ , h = (n -1) p + 1. To minimize the finite sample bias, here, the inverse function of Qn is deduced as Fn (x) := 1 n

x-X cf X cf +1 -X cf + cf , based on Hyndman and Fan's definition, or Fn (x) := 1 n-1 cf -1 +

x-X cf X cf +1 -X cf , based on the latter definition, where cf = n i=1 1 X i ≤x , 1A is the indicator of event A.

The quantile mean uses the location-scale invariant in a different way, as shown in the subsequent proof.

Theorem A.2. qm d= F (µ)-F (WL k,ϵ,γ ) F (WL k,ϵ,γ )- γ 1+γ ,k,ϵ,γ,WL
is a consistent mean estimator for a location-scale distribution provided that the means are finite and F (µ), F (WL k,ϵ,γ ) and γ 1+γ are all within the range of [γϵ, 1 -ϵ], where µ and WL k,ϵ,γ are location parameters from that location-scale distribution. If WL = SWHLM, qm is also consistent for any symmetric distributions.

Proof. When F (WL k,ϵ,γ ) ≥ γ 1+γ , the solution of F (WL k,ϵ,γ ) -γ 1+γ d + F (WL k,ϵ,γ ) = F (µ) is d = F (µ)-F (WL k,ϵ,γ ) F (WL k,ϵ,γ )- γ 1+γ
.

The d value for the case where

F (WL k,ϵ,γ,n ) < γ
1+γ is the same. The definitions of the location and scale parameters are such that they must satisfy

F (x; λ, µ) = F ( x-µ λ ; 1, 0), then F (WL(k, ϵ, γ); λ, µ) = F ( λWL 0 (k,ϵ,γ)+µ-µ λ ; 1, 0) = F (WL0(k, ϵ, γ); 1, 0)
. It follows that the percentile of any weighted L-statistic is free of λ and µ for a location-scale distribution. Therefore d in qm is also invariably a constant. For the symmetric case,

F (SWHLM k,ϵ ) = F (µ) = F (Q( 1 2 )) = 1
2 is valid for any symmetric distribution with a finite second moment, as the same values correspond to same percentiles. Then,

qm d,k,ϵ,SWHLM = F -1 F (SWHLM k,ϵ ) -1 2 d + F (µ) = F -1 (0 + F (µ)) = µ.
To avoid inconsistency due to post-adjustment, F (µ), F (WL k,ϵ,γ ) and γ 1+γ must reside within the range of [γϵ, 1 -ϵ]. All results are now proven.

The cdf of the Pareto distribution is FP ar

(x) = 1 -xm x α .
So, set the d value in qm with two arbitrary percentiles p1 and p2, dP ar,qm =

1-

xm αxm α-1 α -1- xm xm (1-p 1 ) -1 α α 1- xm xm(1-p 1 ) -1 α α -1- xm xm(1-p 2 ) -1 α α = 1-( α-1 α ) α -p 1 p 1 -p 2
. The d value in qm for the exponential distribution is always identical to dP ar,qm as α → ∞,

since limα→∞ α-1 α α = 1
e and the cdf of the exponential distribution is Fexp

(x) = 1 -e -λ -1 x , then dexp,qm = (1-e -1 )-1-e -ln 1 1-p 1 1-e -ln 1 1-p 1 -1-e -ln 1 1-p 2 = 1-1 e -p 1 p 1 -p 2 .
So, for the Weibull, gamma, Pareto, lognormal and generalized Gaussian distribution, qm

d= Fexp (µ)-Fexp (SWHLM k,ϵ ) Fexp(SWHLM k,ϵ )-1 2 ,k,ϵ,SWHLM
is also consistent for at least one particular case, provided that µ and SWHLM k,ϵ are different location parameters from an exponential distribution and F (µ), F (SWHLM k,ϵ ) and 1 2 are all within the range of

[ϵ, 1 -ϵ]. Also let SWHLM k,ϵ,γ = BM ν=3,ϵ= 1 24 and µ = λ, then d = Fexp(µ)-Fexp(BM ν=3,ϵ= 1 24 ) Fexp(BM ν=3,ϵ= 1 24
)- To study the impact of the choice of WLs in rm and qm, it is constructive to recall that a weighted L-statistic is a combination of order statistics. While using a less-biased weighted L-statistic can generally enhance performance (SI Dataset S1), there is a greater risk of violation in the semiparametric framework. However, the mean-WAϵ,γ-γ-median inequality is robust to slight fluctuations of the QA function of the underlying distribution. Suppose for a right-skewed distribution, the QA function is generally decreasing with respect to ϵ in [0, u], but increasing in [u, 1 1+γ ], since all quantile averages with breakdown points from ϵ to 1 1+γ will be included in the computation of WAϵ,γ, as long as 1 1+γ -u ≪ 1 1+γ -γϵ, and other portions of the QA function satisfy the inequality constraints that define the νth γ-orderliness on which the WAϵ,γ is based, if 0 ≤ γ ≤ 1, the mean-WAϵ,γ-γ-median inequality still holds. This is due to the violation of νth γ-orderliness being bounded, when 0 ≤ γ ≤ 1, as shown in REDS I and therefore cannot be extreme for unimodal distributions with finite second moments. For instance, the SQA function of the Weibull distribution is non-monotonic with respect to ϵ when the shape parameter α > -median inequality retains valid when α ≤ 3.387. Another key factor in determining the risk of violation of orderliness is the skewness of the distribution. In REDS I, it was demonstrated that in a family of distributions differing by a skewness-increasing transformation in van Zwet's sense, the violation of orderliness, if it happens, only occurs as the distribution nears symmetry [START_REF] Wr Van Zwet | Convex Transformations of Random Variables: Nebst Stellingen[END_REF]. When γ = 1, the over-corrections in rm and qm are dependent on the SWAϵ-median difference, which can be a reasonable measure of skewness after standardization [START_REF] Al Bowley | Elements of statistics[END_REF][START_REF] Ra Groeneveld | Measuring skewness and kurtosis[END_REF], implying that the over-correction is often tiny with moderate d. This qualitative analysis suggests the general reliability of rm and qm based on the mean-WAϵ,γ-γ-median inequality, especially for unimodal distributions with finite second moments when 0 ≤ γ ≤ 1. Extending this rationale to other weighted L-statistics is possible, since the γ-U -orderliness can also be bounded with certain assumptions, as discussed previously.

Consider two continuous distributions belonging to the same location-scale family, according to Theorem 3 in REDS II, their corresponding kth central moment kernel distributions only differ in scaling. Although strict complete νth orderliness is difficult to prove, following the same logic as discussed above, even if the inequality may be violated in a small range, the mean-SWAϵ-median inequality remains valid, in most cases, for the central moment kernel distribution. Define the recombined kth central moment as

rkm d,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,n,WHLkm 1 ,WHLkm 2 := RI d,h k =ψ k ,k 1 =k,k 2 =k,k 1 ,k 2 ,ϵ 1 ,ϵ 2 ,γ 1 ,γ 2 ,n,LU 1 =WHLkm 1 ,LU 2 =WHLkm 2 .
Then, assuming finite kth central moment and applying the same logic as in Theorem A.1, rkm

d= µ k -WHLkm 1 k 1 ,ϵ 1 ,γ 1 WHLkm 1k 1 ,ϵ 1 ,γ 1 -WHLkm 2 k 2 ,ϵ 2 ,γ 2 ,k 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,WHLkm
qkm d,k,ϵ,γ,n,WHLkm := QI d,h k =ψ k ,k=k,k,ϵ,γ,n,LU =WHLkm . qkm d= F ψ k (µ k )-F ψ k (WHLkm k,ϵ,γ ) F ψ k (WHLkm k,ϵ,γ )- γ 1+γ ,k,ϵ,γ,WHLkm
is also a consistent kth central moment estimator for a location-scale distribution provided that the kth central moment is finite and

F ψ k (µ k ), F ψ k (WHLkm k,ϵ,γ ) and γ
1+γ are all within the range of [γϵ, 1 -ϵ], where µ k and WHLkm k,ϵ,γ are different kth central moment parameters from that location-scale distribution. According to Theorem 2 in REDS II, if the original distribution is unimodal, the central moment kernel distribution is always a heavy-tailed distribution, as the degree term amplifies its skewness and tailedness. From the better performance of the quantile mean in heavy-tailed distributions, the quantile kth central moments are generally better than the recombined kth central moments regarding asymptotic bias.

Finally, the recombined standardized kth moment is defined to be

rskm ϵ=min (ϵ 1 ,ϵ 2 ),k 1 ,k 2 ,k 3 ,k 4 ,γ 1 ,γ 2 ,γ 3 ,γ 4 ,n,WHLkm 1 ,WHLkm 2 ,WHLvar 1 ,WHLvar 2 := rkm d,k 1 ,k 2 ,ϵ 1 ,γ 1 ,γ 2 ,n,WHLkm 1 ,WHLkm 2 (rvar d,k 3 ,k 4 ,ϵ 2 ,γ 3 ,γ 4 ,n,WHLvar 1 ,WHLvar 2 ) k/2 .
The quantile standardized kth moment is defined similarly,

qskm ϵ=min (ϵ 1 ,ϵ 2 ),k 1 ,k 2 ,γ 1 ,γ 2 ,n,WHLkm,WHLvar := qkm d,k 1 ,ϵ 1 ,γ 1 ,n,WHLkm (qvar d,k 2 ,ϵ 2 ,γ 2 ,n,WHLvar ) k/2 .

B. A shape-scale distribution as the consistent distribution.

In the last section, the parametric robust estimation is limited to a location-scale distribution, with the location parameter often being omitted for simplicity. For improved fit to observed skewness or kurtosis, shape-scale distributions with shape parameter (α) and scale parameter (λ) are commonly utilized. Weibull, gamma, Pareto, lognormal, and generalized Gaussian distributions (when µ is a constant) are all shapescale unimodal distributions. Furthermore, if either the shape parameter α or the skewness or kurtosis is constant, the shapescale distribution is reduced to a location-scale distribution.

Let D(|skewness|, kurtosis, k, etype, dtype, n) = d ikm denote the function to specify d values, where the first input is the absolute value of the skewness, the second input is the kurtosis, the third is the order of the central moment (if k = 1, the mean), the fourth is the type of estimator, the fifth is the type of consistent distribution, and the sixth input is the sample size. For simplicity, the last three inputs will be omitted in the following discussion. Hold in awareness that since skewness and kurtosis are interrelated, specifying d values for a shapescale distribution only requires either skewness or kurtosis, while the other may be also omitted. Since many common shape-scale distributions are always right-skewed (if not, only the right-skewed or left-skewed part is used for calibration, while the other part is omitted), the absolute value of the skewness should be the same as the skewness of these distributions. This setting also handles the left-skew scenario well.

For recombined moments up to the fourth ordinal, the object of using a shape-scale distribution as the consistent distribution is to find solutions for the system of equations rvar(WHLvar,γmvar,D(rkurt,2)) 2 . Clearly, this is an overdetermined nonlinear system of equations, given that the skewness and kurtosis are interrelated for a shape-scale distribution. Since an overdetermined system constructed with random coefficients is almost always inconsistent, it is natural to optimize them separately using the fixed-point iteration (see Algorithm 1, only rkurt is provided, others are the same).

                   rm (WHLM, γm, D(|rskew|, rkurt, 1)) = µ rvar (WHLvar, γmvar, D(|rskew|, rkurt, 2)) = µ2 rtm (WHLtm, γmtm, D(|rskew|, rkurt, 3)) = µ3 rf m (WHLf m, γmf m, D(|rskew|, rkurt, 4) = µ4 rskew = µ 3 µ 3 2 2 rkurt = µ 4
i = i + 1 rkurti-1 ← rkurti 6: rkurti ← κ(rkurti-1) until i > maxit or |rkurti -rkurti-1| < δ
▷ maxit is the maximum number of iterations, δ is a small positive number.

The following theorem shows the validity of Algorithm 1. Then, assuming the mkms are all equal to zero and

varD = 1, κ(rkurt D ) rkurt D = f m D -SWHLf m D SWHLf m D (SWHLf m)+SWHLf m rkurt D SWHLvar SWHLvar D 2 = f m D -SWHLf m D SWHLf m D +1 (SWHLf m) f m D SWHLvar SWHLvar D 2 = SWHLf mSWHLvar 2 D SWHLf m D SWHLvar 2 = SWHLf m SWHLvar 2 SWHLf m D SWHLvar D 2 = SWHLkurt SWHLkurt D .
Since SWHLf mD are from the same fourth central moment kernel distribution as f mD = rkurtDvarD 2 , according to the definition of U -congruence, an increase in f mD will also result in an increase in SWHLf mD. Combining with Theorem 3 in REDS II, SWHLkurt is a measure of kurtosis that is invariant to location and scale, so lim rkurt D →∞ κ(rkurt D ) rkurt D < 1. As a result, if there is at least one fixed point, let the largest one be f ixmax, then it is attracting since

| ∂(κ(rkurt D )) ∂(rkurt D ) | < 1 for all rkurtD ∈ [f ixmax, kurtosismax],
where kurtosismax is the maximum kurtosis available in D.

As a result of Theorem B.1, assuming continuity, mkms are all equal to zero, and U -congruence of the central moment kernel distributions, Algorithm 1 converges surely provided that a fixed point exists within the domain of D. At this stage, D can only be approximated through a Monte Carlo study. The continuity of D can be ensured by using linear interpolation. One common encountered problem is that the domain of D depends on both the consistent distribution and the Monte Carlo study, so the iteration may halt at the boundary if the fixed point is not within the domain. However, by setting a proper maximum number of iterations, the algorithm can return the optimal boundary value. For quantile moments, the logic is similar, if the percentiles do not exceed the breakdown point. If this is the case, consistent estimation is impossible, and the algorithm will stop due to the maximum number of iterations. The fixed point iteration is, in principle, similar to the iterative reweighing in Huber M -estimator, but an advantage of this algorithm is that it is solely related to the inputs in Algorithm 1 and is independent of the sample size. Since they are consistent for a shape-scale distribution, |rskew| can specify drm and dtm, rkurt can specify drvar and d rf m . Algorithm 1 enables the robust estimations of all four moments to reach a near-consistent level for common unimodal distributions (Table 1, SI Dataset S1), just using the Weibull distribution as the consistent distribution.

C. Critical points, lines, and Multiple consistent distributions.

In 1895, Pearson considered how to construct probability distributions in which the skewness and kurtosis could be adjusted equally freely. The Pearson distribution is a family of unimodal continuous probability distribution functions that satisfy the following differential equation dP (x) dx = -(a 0 +a 1 x+a 2 x 2 )P (x) b 0 +b 1 x+b 2 x 2 +b 3 x 3 , where P (x) is the pdf of the Pearson distribution, and a0, a1, a2, b0, b1, b2, and b3 are constants that determine the specific type of Pearson distribution. This differential equation ensures that the distribution is unimodal and that the density function is continuous. The Pearson family subsumes many common unimodal distributions, e.g., beta distribution, Cauchy distribution, Chi-squared distribution, continuous uniform distribution, gamma distribution, inverse-chi-squared distribution, inverse-gamma distribution, normal distribution, and tudent's t-distribution. Plotting the kurtosis-skewness lines of the five shape-scale unimodal distributions used in this series into the Pearson diagram. One can immediately identify two critical points, one is 9-2, which is the exponential distribution or the limiting form of the Pareto distribution, another is the normal distribution. This further explains why using the exponential distribution as the consistent distribution has excellent performance even better than using the Weibull distribution as the consistent distribution, since it lays in this critical point. Then, we further consider the critical lines in the Pearson diagram (Figure 1). Consider if the underlying distribution is not the Weibull distribution, but the gamma distribution, then, using the Weibull distribution as an initial guess, calculating the Euclidean distance of the obtained kurtosis-skewness to the kurtosis-skewness lines of other consistent distributions, choosing the one that has the smallest distance, then recalculate the kurtosis-skewness according to the new consistent distribution, iterately, we can use multiple distributions as consistent distributions.

D. Variance.

As one of the fundamental theorems in statistics, the central limit theorem declares that the standard deviation of the limiting form of the sampling distribution of the sample mean is σ √ n . The principle, asymptotic normality, was later applied to the sampling distributions of robust location estimators [START_REF] Newcomb | A generalized theory of the combination of observations so as to obtain the best result[END_REF][START_REF] Daniell | Observations weighted according to order[END_REF][START_REF] Mosteller | On some useful" inefficient" statistics[END_REF][START_REF] Rao | Advanced statistical methods in biometric research[END_REF][START_REF] Bickel | Some contributions to the theory of order statistics[END_REF][START_REF] Chernoff | Asymptotic distribution of linear combinations of functions of order statistics with applications to estimation[END_REF][START_REF] Lecam | On the assumptions used to prove asymptotic normality of maximum likelihood estimates[END_REF][START_REF] Bickel | Descriptive statistics for nonparametric models i[END_REF][START_REF] Bickel | Descriptive statistics for nonparametric models ii[END_REF][START_REF] Janssen | Asymptotic normality for a general class of statistical functions and applications to measures of spread[END_REF]. Daniell (1920) stated ( 13) that comparing the efficiencies of various kinds of estimators is useless unless they all tend to coincide asymptotically. Bickel and Lehmann, also in the landmark series [START_REF] Bickel | Descriptive statistics for nonparametric models i[END_REF][START_REF] Bickel | Descriptive statistics for nonparametric models ii[END_REF], argued that meaningful comparisons of the efficiencies of various kinds of location estimators can be accomplished by studying their standardized variances, asymptotic variances, and efficiency bounds. Standardized variance,

Var( θ) θ 2
, allows the use of simulation studies or empirical data to compare the variances of estimators of distinct parameters. However, a limitation of this approach is the inverse square dependence of the standardized variance on θ. If Var θ1 = Var θ2 , but θ1 is close to zero and θ2 is relatively large, their standardized variances will still differ dramatically. Here, the scaled standard error (SSE) is proposed as a method for estimating the variances of estimators measuring the same attribute, offering a standard error more comparable to that of the sample mean and much less influenced by the magnitude of θ. Definition D.1 (Scaled standard error). Let Ms i s j ∈ R i×j denote the sample-by-statistics matrix, i.e., the first column corresponds to θ, which is the mean or a U -central moment measuring the same attribute of the distribution as the other columns, the second to the jth column correspond to j -1 statistics required to scale, θr 1 , θr 2 , . . ., θr j-1 . Then, the

scaling factor S = 1, θ r 1 θ m , θ r 2 θ m , . . . , θr j-1 θ m T is a j × 1 matrix,
which θ is the mean of the column of Ms i s j . The normalized matrix is M N s i s j = Ms i s j S. The SSEs are the unbiased standard deviations of the corresponding columns of M N s i s j .

The U -central moment (the central moment estimated by using U -statistics) is essentially the mean of the central moment kernel distribution, so its standard error should be generally close to σ km √ n , although not exactly since the kernel distribution is not i.i.d., where σ km is the asymptotic standard deviation of the central moment kernel distribution. If the statistics of interest coincide asymptotically, then the standard errors should still be used, e.g, for symmetric location estimators and odd ordinal central moments for the symmetric distributions, since the scaled standard error will be too sensitive to small changes when they are zero.

The SSEs of all robust estimators proposed here are often, although many exceptions exist, between those of the sample median and those of the sample mean or median central moments and U -central moments (SI Dataset S1). This is because similar monotonic relations between breakdown point and variance are also very common, e.g., Bickel and Lehmann (20) proved that a lower bound for the efficiency of TMϵ to sample mean is (1 -2ϵ) 2 and this monotonic bound holds true for any distribution. However, the direction of monotonicity differs for distributions with different kurtosis. Lehmann andScheffé (1950, 1955) [START_REF] El Lehmann | Completeness, similar regions, and unbiased estimation-part i in Selected works of EL Lehmann[END_REF][START_REF] El Lehmann | Completeness, similar regions, and unbiased estimation-part II[END_REF] in their two early papers provided a way to construct a uniformly minimum-variance unbiased estimator (UMVUE). From that, the sample mean and unbiased sample second moment can be proven as the UMVUEs for the population mean and population second moment for the Gaussian distribution. While their performance for sub-Gaussian distributions is generally satisfied, they perform poorly when the distribution has a heavy tail and completely fail for distributions with infinite second moments. For sub-Gaussian distributions, the variance of a robust location estimator is generally monotonic increasing as its robustness increases, but for heavy-tailed distributions, the relation is reversed. So, unlike bias, the variance-optimal choice can be very different for distributions with different kurtosis.

Due to combinatorial explosion, the bootstrap [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF], introduced by Efron in 1979, is indispensable for computing central moments in practice. In 1981, Bickel and Freedman [START_REF] Bickel | Some asymptotic theory for the bootstrap[END_REF] showed that the bootstrap is asymptotically valid to approximate the original distribution in a wide range of situations, including U -statistics. The limit laws of bootstrapped trimmed U -statistics were proven by Helmers, Janssen, and Veraverbeke (1990) [START_REF] Helmers | Statistics, and System Theory[END_REF]. In REDS I, the advantages of quasi-bootstrap were discussed [START_REF] Rd Richtmyer | A non-random sampling method, based on congruences, for" monte carlo" problems[END_REF][START_REF] Im Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF][START_REF] Ka Do | Quasi-random resampling for the bootstrap[END_REF]. By using quasi-sampling, the impact of the number of repetitions of the bootstrap, or bootstrap size, on variance is very small (SI Dataset S1). An estimator based on the quasi-bootstrap approach can be seen as a complex deterministic estimator that is not only computationally efficient but also statistical efficient. The only drawback of quasi-bootstrap compared to non-bootstrap is that a small bootstrap size can produce additional finite sample bias (SI Text).

E. Root mean square error . [START_REF] Lai | Adaptive choice of mean or median in estimating the center of a symmetric distribution[END_REF] proposed an estimator that adaptively chooses the mean or median in a symmetric distribution and showed that the choice is typically as good as the better of the sample mean and median regarding variance [START_REF] Lai | Adaptive choice of mean or median in estimating the center of a symmetric distribution[END_REF]. Another approach which can be dated back to Laplace (1812) (31) is using wx + (1 -w)mn as a location estimator and w is deduced to achieve optimal variance. Inspired by Lai et al's approach [START_REF] Lai | Adaptive choice of mean or median in estimating the center of a symmetric distribution[END_REF], in this study, for rkurt, there are 364 combinations based on 14 SWHLf ms and 26 SWHLvars (SI Text). Each combination has a root mean square error (RMSE) for a single-parameter distribution, which can be inferred using a Monte Carlo study. For qkurt, there are another 364 combinations, but if the percentiles of quantile moments exceed the breakdown point, that combination is excluded. Then, the combination with the smallest RMSE, calibrated by a two-parameter distribution, is chosen. Similar to Subsection B, let I(kurtosis, dtype, n) = ikurt WHLf m,WHLvar represent these relationships. In this article, the breakdown points of the SWHLMs in SWHLkm were adjusted to ensure the overall breakdown points were 1 24 , as detailed in Theorem F.1). There are two approaches to determine ikurt. The first one is computing all 364+364 rkurt and qkurt, and then, since lim ikurt→∞ I(ikurt) ikurt < 1, the same fix point iteration algorithm as Algorithm 1 can be used to choose the RMSE-optimum combination. The only difference is that unlike D, I is defined to be discontinuous but linear interpolation can also ensure continuity. The second approach is shown in SI Algorithm 2. The RMSEs of these ikurt from the two approaches can be further determined by a Monte Carlo study. Algorithm 1 can also be used to determine the optimum choice among the two approaches. The 364+364 rkurt and qkurt can form a vector, V kurt, where the

Q V kurt ( 1 5 ) to Q V kurt ( 4 5
) can be used to determine the d values of rkms and qkms. The RMSEs of those rkms and qkms can also be estimated by a Monte Carlo study and the estimator with the smallest RMSE of each ordinal is named as ikm. When k is even, the ikurt determined by Ism (detailed in the SI Text) is used to determine ikm. This approach yields results that are often nearly optimal (SI Dateset S1). The estimations of skewness and ikm, when k is odd, follow the same logic.

In general, the variances of invariant central moments are much smaller than those of corresponding unbiased sample central moments (deduced by Cramér [START_REF] Cramér | Mathematical methods of statistics[END_REF][START_REF] Gerlovina | Computer algebra and algorithms for unbiased moment estimation of arbitrary order[END_REF]), except that of the corresponding second central moment (Table 1). F. Robustness. The measure of robustness to gross errors used in this series is the breakdown point proposed by Hampel [START_REF] Hampel | Contributions to the theory of robust estimation[END_REF] in 1968. In REDS I, it has shown that the median of means (MoM) is asymptotically equivalent to the median Hodge-Lehmann mean. Therefore it is also biased for any asymmetric distribution. However, the concentration bound of MoM depends on 1 n [START_REF] Devroye | Sub-gaussian mean estimators[END_REF], it is quite natural to deduce that it is a consistent robust estimator. The concept, sampledependent breakdown point, is defined to avoid ambiguity. Definition F.1 (Sample-dependent breakdown point). The breakdown point of an estimator θ is called sample-dependent if and only if the upper and lower asymptotic breakdown points, which are the upper and lower breakdown points when n → ∞, are zero and the empirical influence function of θ is bounded. For a full formal definition of the empirical influence function, the reader is referred to Devlin, Gnanadesikan and Kettenring (1975)'s paper [START_REF] Sj Devlin | Robust estimation and outlier detection with correlation coefficients[END_REF]. . The second and third tables present the use of the Weibull distribution as the consistent distribution not plus/plus using the lognormal distribution for the odd ordinal moments optimization and the generalized Gaussian distribution for the even ordinal moments optimization. SQM is the robust mean estimator used in recombined/quantile moments. Unbiased sample central moments (var, tm, f m), U -central moments with quasi-bootstrap (var bs , tm bs , f m bs ), and other estimators were used as comparisons. The generalized Gaussian distribution was excluded for He and Fung M -Estimator and Marks percentile estimator, since the logarithmic function does not produce results for negative inputs. The breakdown points of estimators in the second and third table, besides M -estimators and percentile estimator, are all 1 24 . The tables include the average standardized asymptotic bias (ASAB, as n → ∞), root mean square error (RMSE, at n = 5184), average standardized bias (ASB, at n = 5184) and variance (SE ∨ SSE, at n = 5184) of these estimators, all reported in the units of the standard deviations of the distribution or corresponding kernel distributions. W means that the results were weighted by the number of Google Scholar search results on May 30, 2022 (including synonyms). The calibrations of d values and the computations of ASAB, ASB, and SSE were described in Subsection E, F and SI Methods. Detailed results and related codes are available in SI Dataset S1 and GitHub.

Bear in mind that it differs from the "infinitesimal robustness" defined by Hampel, which is related to whether the asymptotic influence function is bounded [START_REF] Hampel | A general qualitative definition of robustness[END_REF][START_REF] Hampel | The influence curve and its role in robust estimation[END_REF][START_REF] Pj Rousseeuw | Robust statistics: the approach based on influence functions[END_REF]. The proof of the consistency of MoM assumes that it is an estimator with a sample-dependent breakdown point since its breakdown point is b 2n , where b is the number of blocks, then limn→∞ b 2n = 0, if b is a constant and any changes in any one of the points of the sample cannot break down this estimator.

For the LU -statistics, the asymptotic upper breakdown points are suggested by the following theorem, which extends the method in Donoho and Huber (1983)'s proof of the breakdown point of the Hodges-Lehmann estimator [START_REF] Dl Donoho | The notion of breakdown point. A festschrift for Erich L. Lehmann[END_REF]. The proof is given in the SI Text.

Theorem F.1. Given a U-statistic associated with a symmetric kernel of degree k. Then, assuming that as n → ∞, k is a constant, the upper breakdown point of the LU -statistic is

1 -(1 -ϵ0) 1 k
, where ϵ0 is the upper breakdown point of the corresponding LL-statistic.

Remark. If k = 1, 1 -(1 -ϵ0) 1 k = ϵ0,
so this formula also holds for the LL-statistic itself. Here, to ensure the breakdown points of all four moments are the same, 1 24 , since ϵ0 = 1 -(1 -ϵ) k , the breakdown points of all LU -statistics for the second, third, and fourth central moment estimations are adjusted as ϵ0 = 47 576 , 1657 13824 , 51935 331776 , respectively. Every statistic is based on certain assumptions. For instance, the sample mean assumes that the second moment of the underlying distribution is finite. If this assumption is violated, the variance of the sample mean becomes infinitely large, even if the population mean is finite. As a result, the sample mean not only has zero robustness to gross errors, but also has zero robustness to departures. To meaningfully compare the performance of estimators under departures from assumptions, it is necessary to impose constraints on these departures. Bound analysis [START_REF] Gauss | Theoria combinationis observationum erroribus minimis obnoxiae[END_REF] is the first approach to study the robustness to departures, i.e., although all estimators can be biased under departures from the corresponding assumptions, but their standardized maximum deviations can differ substantially [START_REF] Devroye | Sub-gaussian mean estimators[END_REF][START_REF] Bieniek | Comparison of the bias of trimmed and winsorized means[END_REF][START_REF] Danielak | Theory & methods: Exact bounds for the bias of trimmed means[END_REF][START_REF] Li | Worst-case range value-at-risk with partial information[END_REF][START_REF] Bernard | Range value-at-risk bounds for unimodal distributions under partial information[END_REF][START_REF] Mathieu | Concentration study of m-estimators using the influence function[END_REF]. In REDS I, it is shown that another way to qualitatively compare the estimators' robustness to departures from the symmetry assumption is constructing and comparing corresponding semiparametric models. While such comparison is limited to a semiparametric model and is not universal, it is still valid for a wide range of parametric distributions. Bound analysis is a more universal approach since they can be deduced by just assuming regularity conditions [START_REF] Devroye | Sub-gaussian mean estimators[END_REF][START_REF] Bieniek | Comparison of the bias of trimmed and winsorized means[END_REF][START_REF] Danielak | Theory & methods: Exact bounds for the bias of trimmed means[END_REF][START_REF] Li | Worst-case range value-at-risk with partial information[END_REF][START_REF] Mathieu | Concentration study of m-estimators using the influence function[END_REF]. However, bounds are often hard to deduce for complex estimators. Also, sometimes there are discrepancies between maximum bias and average bias. Since the estimators proposed here are all consistent under certain assumptions, measuring their biases is also a convenient way of measuring the robustness to departures. Average standardized asymptotic bias is thus defined as follows. estimator θ is not classified as an RI-statistic, QI-statistic, or LU -statistic, the corresponding U -statistic, which measures the same attribute of the distribution, is utilized to determine the value of σ. For a two-parameter distribution, the first step is setting the lower bound of the kurtosis range of interest μ4 l , the spacing δ, and the bin count C. Then, the average standardized asymptotic bias is defined as

ASAB θ := 1 C δ+ μ4 l ≤ μ4 ≤Cδ+ μ4 l μ4 is a multiple of δ E θ| μ4 θ -θ σ
where μ4 is the kurtosis specifying the two-parameter distribution, E θ| μ4 denotes the expected value given fixed μ4.

Standardization plays a crucial role in comparing the performance of estimators across different distributions. Currently, several options are available, such as using the root mean square deviation from the mode (as in Gauss (1)), the mean absolute deviation, or the standard deviation. The standard deviation is preferred due to its central role in standard error estimation. In Table 1, δ = 0.1, C = 70. For the Weibull, gamma, lognormal and generalized Gaussian distributions, μ4 l = 3 (there are two shape parameter solutions for the Weibull distribution, the lower one is used here). For the Pareto distribution, μ4 l = 9. To provide a more practical and straightforward illustration, all results from five distributions are further weighted by the number of Google Scholar search results. Within the range of kurtosis setting, nearly all WLs and WHLkms proposed here reach or at least come close to their maximum biases (SI Dataset S1). The pseudomaximum bias is thus defined as the maximum value of the biases within the range of kurtosis setting for all five unimodal distributions. In most cases, the pseudo-maximum biases of invariant moments occur in lognormal or generalized Gaussian distributions (SI Dataset S1), since besides unimodality, the Weibull distribution differs entirely from them. Interestingly, the asymptotic biases of TM ϵ= 1 8 and WM ϵ= 1 8

, after averaging and weighting, are 0.107σ and 0.066σ, respectively, in line with the sharp bias bounds of TM2,14:15 and WM2,14:15 (a different subscript is used to indicate a sample size of 15, with the removal of the first and last order statistics), 0.173σ and 0.126σ, for any distributions with finite moments [START_REF] Bieniek | Comparison of the bias of trimmed and winsorized means[END_REF][START_REF] Danielak | Theory & methods: Exact bounds for the bias of trimmed means[END_REF].

Discussion

Statistics, encompassing the collection, analysis, interpretation, and presentation of data, has evolved over time, with various approaches emerging to meet challenges in practice. Among these approaches, the use of probability models and measures of random variables for data analysis is often considered the core of statistics. While the early development of statistics was focused on parametric methods, there were two main approaches to point estimation. The Gauss-Markov theorem [START_REF] Gauss | Theoria combinationis observationum erroribus minimis obnoxiae[END_REF]46) states the principle of minimum variance unbiased estimation which was enriched by [START_REF] Neyman | On the two different aspects of the representative method: The method of stratified sampling and the method of purposive selection[END_REF] [START_REF] Neyman | On the two different aspects of the representative method: The method of stratified sampling and the method of purposive selection[END_REF], [START_REF] Rao | Information and accuracy attainable in the estimation of statistical parameters[END_REF] [START_REF] Rao | Information and accuracy attainable in the estimation of statistical parameters[END_REF], [START_REF] Blackwell | Conditional expectation and unbiased sequential estimation[END_REF] [START_REF] Blackwell | Conditional expectation and unbiased sequential estimation[END_REF], and Lehmann and Scheffé (1950,1955) [START_REF] El Lehmann | Completeness, similar regions, and unbiased estimation-part i in Selected works of EL Lehmann[END_REF][START_REF] El Lehmann | Completeness, similar regions, and unbiased estimation-part II[END_REF]. Maximum likelihood was first introduced by Fisher in 1922 [START_REF] Fisher | On the mathematical foundations of theoretical statistics[END_REF] in a multinomial model and later generalized by Cramér (1946), Hájek (1970), and Le Cam (1972) [START_REF] Lecam | On the assumptions used to prove asymptotic normality of maximum likelihood estimates[END_REF][START_REF] Cramér | Mathematical methods of statistics[END_REF][START_REF] Hájek | Local asymptotic minimax and admissibility in estimation[END_REF]. In 1939, Wald (52) combined these two principles and suggested the use of minimax estimates, which involve choosing an estimator that minimizes the maximum possible loss. Following Huber's seminal work [START_REF] Huber | Robust estimation of a location parameter[END_REF], M -statistics have dominated the field of parametric robust statistics for over half a century. Nonparametric methods, e.g., the Kolmogorov-Smirnov test, Mann-Whitney-Wilcoxon Test, and Hoeffding's independence test, emerged as popular alternatives to parametric methods in 1950s, as they do not make specific assumptions about the underlying distribution of the data. In 1963, Hodges and Lehmann proposed a class of robust location estimators based on the confidence bounds of rank tests [START_REF] Hodges | Estimates of location based on rank tests[END_REF]. In REDS I, when compared to other semiparametric mean estimators with the same breakdown point, the H-L estimator was shown to be the bias-optimal choice, which aligns Devroye, and Lerasle, Lugosi, and Oliveira's conclusion that the median of means is near-optimal in terms of concentration bounds [START_REF] Devroye | Sub-gaussian mean estimators[END_REF] as discussed. The formal study of semiparametric models was initiated by Stein (54) in 1956. Bickel, in 1982, simplified the general heuristic necessary condition proposed by Stein [START_REF] Stein | Efficient nonparametric testing and estimation in Proceedings of the third Berkeley symposium on mathematical statistics and probability[END_REF] and derived sufficient conditions for this type of problem, adaptive estimation [START_REF] Bickel | On adaptive estimation[END_REF]. These conditions were subsequently applied to the construction of adaptive estimates [START_REF] Bickel | On adaptive estimation[END_REF]. It has become increasingly apparent that, in robust statistics, many estimators previously called "nonparametric" are essentially semiparametric as they are partly, though not fully, characterized by some interpretable Euclidean parameters. This approach is particularly useful in situations where the data do not conform to a simple parametric distribution but still have some structure that can be exploited. In 1984, Bickel addressed the challenge of robustly estimating the parameters of a linear model while acknowledging the possibility that the model may be invalid but still within the confines of a larger model [START_REF] Bickel | Parametric robustness: small biases can be worthwhile[END_REF]. He showed by carefully designing the estimators, the biases can be very small. The paradigm shift here opens up the possibility that by defining a large semiparametric model and constructing estimators simultaneously for two or more very different semiparametric/parametric models within the large semiparametric model, then even for a parametric model belongs to the large semiparametric model but not to the semiparametric/parametric models used for calibration, the performance of these estimators might still be near-optimal due to the common nature shared by the models used by the estimators. Maybe it can be named as comparametrics. Closely related topics are "mixture model" and "constraint defined model," which were generalized in Bickel, Klaassen, Ritov, and Wellner's classic semiparametric textbook (1993) [START_REF] Bickel | Efficient and adaptive estimation for semiparametric models[END_REF] and the method of sieves, introduced by Grenander in 1981 [START_REF] Grenander | Abstract Inference[END_REF]. As the building blocks of statistics, invariant moments can reduce the overall errors of statistical results across studies and thus can enhance the replicability of the whole community [START_REF] Jt Leek | Reproducible research can still be wrong: adopting a prevention approach[END_REF][START_REF]Reproducibility and Replicability in Science[END_REF].

Methods

Methods of generating the Table 1 are summarized below, with details in the SI Text. The d values for the invariant moments of the Weibull distribution were approximated using a Monte Carlo study, with the formulae presented in Theorem A.1 and A.2. The computation of I functions is summarized in Subsection E and further explained in the SI Text. The computation of ASABs and ASBs is described in Subsection F. The SEs and SSEs were computed by approximating the sampling distribution using 1000 pseudorandom samples for n = 5184 and 50 pseudorandom samples for n = 2654208. The impact of the bootstrap size, ranging from n = 2.7 × 10 2 to n = 2.765 × 10 4 , on the variance of invariant moments and U -central moments was studied using the SEs and SSEs methods described above. A brute force approach was used to estimate the maximum biases of the robust estimators discussed for the five unimodal distributions. The validity of this approach is discussed in the SI Text.

1 1 -

 1 ln(2) ≈ 3.259 as shown in the SI Text of REDS I, the violation of the second and third orderliness starts near this parameter as well, yet the mean-BM ν=3,ϵ= 1 24

Li PNAS |

  February 17, 2024 | vol. XXX | no. XX | 3

µ 2 2 ,

 2 where µ2, µ3 and µ4 are the population second, third and fourth central moments.|rskew| and rkurt should be the invariant points of the functions ς(|rskew|) = rtm(WHLtm,γmtm,D(|rskew|,3)) rvar(WHLvar,γmvar,D(|rskew|,2))

  ) = rf m(WHLf m,γmf m,D(rkurt,4))

Algorithm 1

 1 rkurt for a shape-scale distribution Input: D; WHLvar; WHLf m; γmvar; γmf m; maxit; δ Output: rkurti-1 i = 0 2: rkurti ← κ(kurtosismax) ▷ Using the maximum kurtosis available in D as an initial guess. repeat 4:

Theorem B. 1 .

 1 Assuming γ = 1 and mkms, where 2 ≤ k ≤ 4, are all equal to zero, |rskew| and rkurt, defined as the largest attracting fixed points of the functions ς(|rskew|) and κ(rkurt), are consistent estimators of μ3 and μ4 for a shape-scale distribution whose kth central moment kernel distributions are U -congruent, as long as they are within the domain of D, where μ3 and μ4 are the population skewness and kurtosis, respectively. Proof. Without loss of generality, only rkurt is considered, while the logic for |rskew| is the same. Additionally, the second central moments of the underlying sample distribution and consistent distribution are assumed to be 1, with other cases simply multiplying a constant factor according to Theorem 3 in REDS II. From the definition of D, κ(rkurt D ) rkurt D = f m D -SWHLf m D SWHLf m D -mf m D (SWHLf m-mf m)+SWHLf m rkurt D var D -SWHLvar D SWHLvar D -mvar D (SWHLvar-mvar)+SWHLvar 2 , where the subscript D indicates that the estimates are from the central moment kernel distributions generated from the consistent distribution, while other estimates are from the underlying distribution of the sample.

Li PNAS | February 17 ,

 17 2024 | vol. XXX | no. XX | 5

Definition F. 2 (

 2 Average standardized asymptotic bias). For a single-parameter distribution, the average standardized asymptotic bias (ASAB) is given by | θ-θ| σ , where θ represents the estimation of θ, and σ denotes the standard deviation of the kernel distribution associated with the LU -statistic. If the Li PNAS | February 17, 2024 | vol. XXX | no. XX | 7

  only the asymptotic behavior is considered, in the absence of subscripts, no assumptions are made. The subsequent theorem shows the significance of a recombined I-statistic.

	Theorem	A.1.	Define	the	recombined	mean
	as

  ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,WL 1 ,WL 2

	Proof. Finding	d	that	make
	rm d,k 1 a	consistent
	mean estimator is equivalent to finding the so-
	lution	of	rm d,k

1 ,WL 2 is a consistent mean estimator for a location-scale distribution, where µ, WL1 k 1 ,ϵ 1 ,γ 1 , and WL2 k 2 ,ϵ 2 ,γ 2 are different location parameters from that location-scale distribution. If γ1 = γ2 = 1, WL = SWHLM, rm is also consistent for any symmetric distributions. 1 ,k 2 ,ϵ=min (ϵ 1 ,ϵ 2 ),γ 1 ,γ 2 ,WL 1 ,WL 2 = µ.

  both consistent mean estimators for any symmetric distribution and the exponential distribution with finite second moments. It's obvious that the asymptotic breakdown points of rm d≈0.103,ν=3,ϵ= 1

				1	=
	-e -1 +e -1+ln 1 -1+ln 2 -e 101898752449325 √ 26068394603446272 6 26068394603446272 6 √ 391 5/6 101898752449325 7 247 √ 3 √ 5 11 √ 7 3 √ 11 247 391 5/6 101898752449325 √ 5 5 6 √ 247 7 391 5/6 26068394603446272 3 √ 11e -1 e 1 2 -101898752449325 √ 5 6 √ 247 26068394603446272 3 11e √ 7 391 5/6 ≈	0.088.	2 Fexp(µ), =
	Fexp(BM ν=3,ϵ= 1 24 [ 1 24 , 23 24 ]. qm d≈0.088,ν=3,ϵ= 1 ) and 1 2 are all within the range of 24 ,BM works better in the fat-tail
	scenarios (SI Dataset S1). Theorem A.1 and A.2 show
	that rm d≈0.103,ν=3,ϵ= 1 24 ,BM,m	and qm d≈0.088,ν=3,ϵ= 1 24 ,BM

are 24 ,BM,m and qm d≈0.088,ν=3,ϵ= 1 24 ,BM are both

1 24

. Therefore they are all invariant means.

  1 ,WHLkm 2 is a consistent kth central moment estimator for a location-scale distribution, where µ k , WHLkm1 k 1 ,ϵ 1 ,γ 1 , and WHLkm2 k 2 ,ϵ 2 ,γ 2 are different kth central moment parameters from that location-scale distribution. Similarly, the quantile will not change after scaling. The quantile kth central moment is thus defined as

Table 1 . Evaluation of invariant moments for five common unimodal distributions in comparison with current popular methods

 1 The first table presents the use of the exponential distribution as the consistent distribution for five common unimodal distributions: Weibull, gamma, Pareto, lognormal, and generalized Gaussian distributions. Popular robust mean estimators discussed in REDS 1 were used as comparisons. The breakdown points of mean estimators in the first table, besides H-L estimator and Huber M -estimator, are all1 8

	Errors	x	TM	H-L	SM	HM	WM	SQM	BM	MoM	MoRM	mHLM	rmexp,BM	qmexp,BM
	WASAB	0.000	0.107	0.088	0.078	0.078	0.066	0.048	0.048	0.034	0.035	0.034	0.002		0.003
	WRMSE	0.014	0.111	0.092	0.083	0.083	0.070	0.053	0.053	0.041	0.041	0.038	0.017		0.018
	WASBn=5184	0.000	0.108	0.089	0.078	0.079	0.066	0.048	0.048	0.034	0.036	0.033	0.002		0.003
	WSE ∨ WSSE	0.014	0.014	0.014	0.015	0.014	0.014	0.014	0.015	0.017	0.014	0.014	0.017		0.017
	Errors		HFMµ	MPµ	rm	qm	im	var	var bs	Tsd 2	HFMµ 2	MPµ 2	rvar	qvar	ivar
	WASAB		0.037	0.043	0.001	0.002	0.001	0.000	0.000	0.200	0.027	0.042	0.005	0.018	0.003
	WRMSE		0.049	0.055	0.015	0.015	0.014	0.017	0.017	0.198	0.042	0.062	0.019	0.026	0.019
	WASBn=5184		0.038	0.043	0.001	0.002	0.001	0.000	0.001	0.198	0.027	0.043	0.005	0.018	0.003
	WSE ∨ WSSE		0.018	0.021	0.015	0.015	0.014	0.017	0.017	0.015	0.024	0.032	0.018	0.017	0.018
	Errors	tm	tm bs	HFMµ 3	MPµ 3	rtm	qtm	itm	f m	f m bs	HFMµ 4	MPµ 4	rf m	qf m	if m
	WASAB	0.000	0.000	0.052	0.059	0.006	0.083	0.034	0.000	0.000	0.037	0.046	0.024	0.038	0.011
	WRMSE	0.019	0.018	0.063	0.074	0.018	0.083	0.044	0.026	0.023	0.049	0.062	0.037	0.043	0.029
	WASBn=5184	0.001	0.003	0.052	0.059	0.007	0.082	0.038	0.001	0.009	0.037	0.047	0.024	0.036	0.013
	WSE ∨ WSSE	0.019	0.018	0.021	0.091	0.015	0.012	0.017	0.024	0.021	0.020	0.027	0.021	0.020	0.022

To whom correspondence should be addressed. E-mail: tl@biomathematics.org PNAS | February 17,

| vol. XXX | no. XX | 1-9

| Li

| Li

| Li

Data and Software Availability. Data for Table 1 are given in SI Dataset S1-S4. All codes have been deposited in GitHub.