Robust estimations from distribution structures: I. Mean

I n 1823, Gauss (1) proved that for any unimodal distribution, |m -µ| ≤ 3 4 ω and σ ≤ ω ≤ 2σ, where µ is the population mean, m is the population median, ω is the root mean square deviation from the mode, and σ is the population standard deviation. This pioneering work revealed that, the potential bias of the median with respect to the mean is bounded in units of a scale parameter under certain assumptions. In 2018, Li, Shao, Wang, and Yang (2) proved the bias bound of any quantile for arbitrary continuous distributions with finite second moments. [START_REF] Bernard | Range value-at-risk bounds for unimodal distributions under partial information[END_REF] (3) further refined these bounds for unimodal distributions with finite second moments and extended to the bounds of symmetric quantile averages. They showed that m has the smallest maximum distance to µ among all symmetric quantile averages (SQA ϵ ). Daniell, in 1920, (4) analyzed a class of estimators, linear combinations of order statistics, and identified that the ϵ-symmetric trimmed mean (STMϵ) belongs to this class. Another popular choice, the ϵ-symmetric Winsorized mean (SWMϵ), named after Winsor and introduced by Tukey (5) and Dixon [START_REF] Wj Dixon | Simplified Estimation from Censored Normal Samples[END_REF] in 1960, is also an L-estimator. [START_REF] Bieniek | Comparison of the bias of trimmed and winsorized means[END_REF] derived exact bias upper bounds of the Winsorized mean based on Danielak and Rychlik's work (2003) on the trimmed mean for any distribution with a finite second moment and confirmed that the former is smaller than the latter [START_REF] Danielak | Theory & methods: Exact bounds for the bias of trimmed means[END_REF][START_REF] Bieniek | Comparison of the bias of trimmed and winsorized means[END_REF]. [START_REF] Oliveira | The sub-gaussian property of trimmed means estimators[END_REF] and [START_REF] Lugosi | Robust multivariate mean estimation: The optimality of trimmed mean[END_REF] [START_REF] Oliveira | The sub-gaussian property of trimmed means estimators[END_REF][START_REF] Lugosi | Robust multivariate mean estimation: The optimality of trimmed mean[END_REF] derived the concentration bound of the trimmed mean. In 1963, Hodges and Lehmann [START_REF] Hodges | Estimates of location based on rank tests[END_REF] proposed a class of nonparametric location estimators based on rank tests and, from the Wilcoxon signed-rank statistic (12), deduced the median of pairwise means as a robust location estimator for a symmetric population. Both L-statistics and R-statistics achieve robustness essentially by removing a certain proportion of extreme values, therefore, they have predefined breakdown points [START_REF] Hampel | Contributions to the theory of robust estimation[END_REF]. In 1964, Huber [START_REF] Huber | Robust estimation of a location parameter[END_REF] generalized maximum likeli-hood estimation to the minimization of the sum of a specific loss function, which measures the residuals between the data points and the model's parameters. Some L-estimators are also M -estimators, e.g., the sample mean is an M -estimator with a squared error loss function, the sample median is an M -estimator with an absolute error loss function [START_REF] Huber | Robust estimation of a location parameter[END_REF]. The Huber M -estimator is obtained by applying the Huber loss function that combines elements of both squared error and absolute error to achieve robustness against gross errors and high efficiency for contaminated Gaussian distributions [START_REF] Huber | Robust estimation of a location parameter[END_REF]. Sun, Zhou, and Fan (2020) examined the concentration bounds of the Huber M -estimator [START_REF] Sun | Adaptive huber regression[END_REF]. In 2012, Catoni proposed an M -estimator for heavy-tailed samples with finite variance [START_REF] Catoni | Challenging the empirical mean and empirical variance: a deviation study[END_REF]. Xu (2021) [START_REF] Chen | A generalized Catoni's M-estimator under finite α-th moment assumption with α ∈ (1, 2)[END_REF] proposed a generalized Catoni M -estimator and showed that it has a better worse-case performance than the empirical mean. [START_REF] Mathieu | Concentration study of m-estimators using the influence function[END_REF] [START_REF] Mathieu | Concentration study of m-estimators using the influence function[END_REF] further derived the concentration bounds of M -estimators and demonstrated that, by selecting the tuning parameter which depends on the variance, these M -estimator can also be a sub-Gaussian estimator. The concept of the median of means (MoM k,b= n k ,n ) was first introduced by Nemirovsky and Yudin (1983) in their work on stochastic optimization [START_REF] As Nemirovskij | Problem complexity and method efficiency in optimization[END_REF], while later was revisited in Jerrum, [START_REF] Mr Jerrum | Random generation of combinatorial structures from a uniform distribution[END_REF], [START_REF] Mr Jerrum | Random generation of combinatorial structures from a uniform distribution[END_REF] and Alon, Matias and Szegedy (1996) [START_REF] Alon | The space complexity of approximating the frequency moments[END_REF]'s works. Given its good performance even for distributions with infinite second moments, the MoM has received increasing attention over the past decade [START_REF] Hsu | Heavy-tailed regression with a generalized median-of-means in International Conference on Machine Learning[END_REF][START_REF] Devroye | Sub-gaussian mean estimators[END_REF][START_REF] Laforgue | On medians of (randomized) pairwise means in International Conference on Machine Learning[END_REF][START_REF] Lecué | Robust machine learning by median-of-means: Theory and practice[END_REF]. Devroye, Lerasle, Lugosi, and [START_REF] Devroye | Sub-gaussian mean estimators[END_REF] showed that MoM k,b= n k ,n nears the optimum of sub-Gaussian mean estimation with regards to concentration bounds when the distribution has a heavy tail [START_REF] Devroye | Sub-gaussian mean estimators[END_REF]. Laforgue, Clemencon, and Bertail (2019) proposed the median of randomized means (MoRM k,b,n ) [START_REF] Laforgue | On medians of (randomized) pairwise means in International Conference on Machine Learning[END_REF], wherein, rather than partitioning, an ar-

Significance Statement

In 1964, van Zwet introduced the convex transformation order for comparing the skewness of two distributions. This paradigm shift played a fundamental role in defining robust measures of distributions, from spread to kurtosis. Here, instead of examining the stochastic ordering between two distributions, the orderliness of quantile combinations within a distribution is investigated. By classifying distributions through the signs of derivatives, two series of sophisticated robust mean estimators are deduced. Nearly all common nonparametric robust location estimators are found to be special cases thereof.

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

The author declares no competing interest.

D R A F T

bitrary number, b, of blocks are built independently from the sample, and showed that MoRM k,b,n has a better nonasymptotic sub-Gaussian property compared to MoM k,b= n k ,n . In fact, asymptotically, the Hodges-Lehmann (H-L) estimator is equivalent to MoM k=2,b= n k and MoRM k=2,b , and they can be seen as the pairwise mean distribution is approximated by the sampling without replacement and bootstrap, respectively. When k ≪ n, the difference between sampling with replacement and without replacement is negligible. For the asymptotic validity, readers are referred to the foundational works of [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF] [START_REF] Efron | Bootstrap methods: Another look at the jackknife[END_REF], [START_REF] Bickel | On adaptive estimation[END_REF]Freedman (1981, 1984) [START_REF] Bickel | Some asymptotic theory for the bootstrap[END_REF][START_REF] Bickel | Asymptotic normality and the bootstrap in stratified sampling[END_REF], and Helmers, Janssen, and Veraverbeke (1990) [START_REF] Helmers | Statistics, and System Theory[END_REF].

Here, the ϵ,b-stratified mean is defined as

SM ϵ,b,n := b n    b-1 2bϵ j=1 (2bj-b+1)nϵ b-1 i j = (2bj-b-1)nϵ b-1 +1 Xi j    ,
where X1 ≤ . . . ≤ Xn denote the order statistics of a sample of n independent and identically distributed random variables X1, . . ., Xn. b ∈ N, b ≥ 3, and b mod 2 = 1. The definition was further refined to guarantee the continuity of the breakdown point by incorporating an additional block in the center when ⌊ b-1 2bϵ ⌋ mod 2 = 0, or by adjusting the central block when

⌊ b-1 2bϵ ⌋ mod 2 = 1 (SI Text). If the subscript n is omitted, only the asymptotic behavior is considered. If b is omitted, b = 3 is assumed. SM ϵ,b=3 is equivalent to STMϵ, when ϵ > 1 6 . When b-1 2ϵ ∈ N, the basic idea of the stratified mean is to distribute the data into b-1
2ϵ equal-sized non-overlapping blocks according to their order. Then, further sequentially group these blocks into b equal-sized strata and compute the mean of the middle stratum, which is the median of means of each stratum. In situations where i mod 1 ̸ = 0, a potential solution is to generate multiple smaller samples that satisfy the equality by sampling without replacement, and subsequently calculate the mean of all estimations. The details of determining the smaller sample size and the number of sampling times are provided in the SI Text. Although the principle resembles that of the median of means, SM ϵ,b,n is different from MoM k= n b ,b,n as it does not include the random shift. Additionally, the stratified mean differs from the mean of the sample obtained through stratified sampling methods, introduced by Neyman (1934) [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] or ranked set sampling [START_REF] Mcintyre | A method for unbiased selective sampling, using ranked sets[END_REF], introduced by McIntyre in 1952, as these sampling methods aim to obtain more representative samples or improve the efficiency of sample estimates, but the sample means based on them are not robust. When b mod 2 = 1, the stratified mean can be regarded as replacing the other equal-sized strata with the middle stratum, which, in principle, is analogous to the Winsorized mean that replaces extreme values with less extreme percentiles. Furthermore, while the bounds confirm that the Winsorized mean and median of means outperform the trimmed mean [START_REF] Danielak | Theory & methods: Exact bounds for the bias of trimmed means[END_REF][START_REF] Bieniek | Comparison of the bias of trimmed and winsorized means[END_REF][START_REF] Devroye | Sub-gaussian mean estimators[END_REF] in worst-case performance, the complexity of bound analysis makes it difficult to achieve a complete and intuitive understanding of these results. Also, a clear explanation for the average performance of them remains elusive. The aim of this paper is to define a series of semiparametric models using the signs of derivatives, reveal their elegant interrelations and connections to parametric models, and show that by exploiting these models, two sets of sophisticated mean estimators can be deduced, which exhibit strong robustness to departures from assumptions.

Quantile Average and Weighted Average

The symmetric trimmed mean, symmetric Winsorized mean, and stratified mean are all L-estimators. More specifically, they are symmetric weighted averages, which are defined as

SWAϵ,n := ⌈ n 2 ⌉ i=1 X i +X n-i+1 2 wi ⌈ n 2 ⌉ i=1 wi
, where wis are the weights applied to the symmetric quantile averages according to the definition of the corresponding Lestimators. For example, for the ϵ-symmetric trimmed mean, wi = 0, i < nϵ 1, i ≥ nϵ , when nϵ ∈ N. The mean and median are indeed two special cases of the symmetric trimmed mean.

To extend the symmetric quantile average to the asymmetric case, two definitions for the ϵ,γ-quantile average (QA ϵ,γ,n ) are proposed. The first definition is:

1 2 ( Qn(γϵ) + Qn(1 -ϵ)), [1]
and the second definition is:

1 2 ( Qn(ϵ) + Qn(1 -γϵ)), [2]
where Qn(p) is the empirical quantile function; γ is used to adjust the degree of asymmetry, γ ≥ 0; and 0 ≤ ϵ ≤ 1 1+γ . For trimming from both sides, [START_REF] Gauss | Theoria combinationis observationum erroribus minimis obnoxiae[END_REF] and [START_REF] Li | Worst-case range value-at-risk with partial information[END_REF] are essentially equivalent. The first definition along with γ ≥ 0 and 0 ≤ ϵ ≤ 1 1+γ are assumed in the rest of this article unless otherwise specified, since many common asymmetric distributions are right-skewed, and [START_REF] Gauss | Theoria combinationis observationum erroribus minimis obnoxiae[END_REF] allows trimming only from the right side by setting γ = 0.

Analogously, the weighted average can be defined as WAϵ,γ,n :=

1 1+γ 0 QA (ϵ0, γ, n) w(ϵ0)dϵ0 1 1+γ 0 w(ϵ0)dϵ0
. For any weighted average, if γ is omitted, it is assumed to be 1. The ϵ, γ-trimmed mean (TMϵ,γ,n) is a weighted average with a left trim size of nγϵ and a right trim size of nϵ, where w(ϵ0) = 0, ϵ0 < ϵ 1, ϵ0 ≥ ϵ . Using this definition, regardless of whether nγϵ / ∈ N or nϵ / ∈ N, the TM computation remains the same, since this definition is based on the empirical quantile function. However, in this article, considering the computational cost in practice, non-asymptotic definitions of various types of weighted averages are primarily based on order statistics. Unless stated otherwise, the solution to their decimal issue is the same as that in SM.

Furthermore, for weighted averages, separating the breakdown point into upper and lower parts is necessary.

Definition .1 (Upper/lower breakdown point). The upper breakdown point is the breakdown point generalized in Davies and Gather (2005)'s paper [START_REF] Davies | Breakdown and groups[END_REF]. The finite-sample upper breakdown point is the finite sample breakdown point defined by Donoho and Huber (1983) [START_REF] Dl Donoho | The notion of breakdown point. A festschrift for Erich L. Lehmann[END_REF] and also detailed in [START_REF] Davies | Breakdown and groups[END_REF]. The (finite-sample) lower breakdown point is replacing the infinity symbol in these definitions with negative infinity.

D R A F T

Classifying Distributions by the Signs of Derivatives

Let P R denote the set of all continuous distributions over R and P X denote the set of all discrete distributions over a countable set X. The default of this article will be on the class of continuous distributions, P R . However, it's worth noting that most discussions and results can be extended to encompass the discrete case, P X , unless explicitly specified otherwise. Besides fully and smoothly parameterizing them by a Euclidean parameter or merely assuming regularity conditions, there exist additional methods for classifying distributions based on their characteristics, such as their skewness, peakedness, modality, and supported interval. In 1956, Stein initiated the study of estimating parameters in the presence of an infinitedimensional nuisance shape parameter [START_REF] Stein | Efficient nonparametric testing and estimation in Proceedings of the third Berkeley symposium on mathematical statistics and probability[END_REF] and proposed a necessary condition for this type of problem, a contribution later explicitly recognized as initiating the field of semiparametric statistics [START_REF] Bickel | On adaptive estimation[END_REF]. In 1982, Bickel simplified Stein's general heuristic necessary condition [START_REF] Stein | Efficient nonparametric testing and estimation in Proceedings of the third Berkeley symposium on mathematical statistics and probability[END_REF], derived sufficient conditions, and used them in formulating adaptive estimates [START_REF] Bickel | On adaptive estimation[END_REF]. A notable example discussed in these groundbreaking works was the adaptive estimation of the center of symmetry for an unknown symmetric distribution, which is a semiparametric model. In 1993, Bickel, Klaassen, Ritov, and Wellner published an influential semiparametrics textbook [START_REF] Bickel | Efficient and adaptive estimation for semiparametric models[END_REF], which categorized most common statistical models as semiparametric models, considering parametric and nonparametric models as two special cases within this classification. Yet, there is another old and commonly encountered class of distributions that receives little attention in semiparametric literature: the unimodal distribution. It is a very unique semiparametric model because its definition is based on the signs of derivatives, i.e., (f ′ (x) > 0 for x ≤ M ) ∧ (f ′ (x) < 0 for x ≥ M ), where f (x) is the probability density function (pdf) of a random variable X, M is the mode. Let PU denote the set of all unimodal distributions. There was a widespread misbelief that the median of an arbitrary unimodal distribution always lies between its mean and mode until Runnenburg (1978) and van [START_REF] Zwet | Mean, median, mode ii[END_REF] [START_REF] Jt Runnenburg | Mean, median, mode[END_REF][START_REF] Zwet | Mean, median, mode ii[END_REF] endeavored to determine sufficient conditions for the mean-median-mode inequality to hold, thereby implying the possibility of its violation. The class of unimodal distributions that satisfy the mean-median-mode inequality constitutes a subclass of PU , denoted by PMMM ⊊ PU . To further investigate the relations of location estimates within a distribution, the γ-orderliness for a right-skewed distribution is defined as

∀0 ≤ ϵ1 ≤ ϵ2 ≤ 1 1 + γ , QA(ϵ1, γ) ≥ QA(ϵ2, γ).
The necessary and sufficient condition below hints at the relation between the mean-median-mode inequality and the γ-orderliness.

Theorem .1. A distribution is γ-ordered if and only if its pdf satisfies the inequality

f (Q(γϵ)) ≥ f (Q(1 -ϵ)) for all 0 ≤ ϵ ≤ 1 1+γ or f (Q(γϵ)) ≤ f (Q(1 -ϵ)) for all 0 ≤ ϵ ≤ 1 1+γ .
Proof. Without loss of generality, consider the case of rightskewed distribution. From the above definition of γ-orderliness, it is deduced that

Q(γϵ-δ)+Q(1-ϵ+δ) 2 ≥ Q(γϵ)+Q(1-ϵ) 2 ⇔ Q(γϵ- δ) -Q(γϵ) ≥ Q(1 -ϵ) -Q(1 -ϵ + δ) ⇔ Q ′ (1 -ϵ) ≥ Q ′ (γϵ),
where δ is an infinitesimal positive quantity. Observing that the quantile function is the inverse function of the cumulative

distribution function (cdf), Q ′ (1 -ϵ) ≥ Q ′ (γϵ) ⇔ F ′ (Q(γϵ)) ≥ F ′ (Q(1-ϵ))
, thereby completing the proof, since the derivative of cdf is pdf.

According to Theorem .1, if a probability distribution is right-skewed and monotonic decreasing, it will always be γordered. For a right-skewed unimodal distribution, if Q(γϵ) > M , then the inequality f (Q(γϵ)) ≥ f (Q(1 -ϵ)) holds. The principle is extendable to unimodal-like distributions. Suppose there is a right-skewed unimodal-like distribution with the first mode, denoted as M1, having the greatest probability density, while there are several smaller modes located towards the higher values of the distribution. Furthermore, assume that this distribution follows the mean-γ-median-first mode inequality, and the γ-median, Q( γ 1+γ ), falling within the first dominant mode (i.e., if

x > Q( γ 1+γ ), f (Q( γ 1+γ )) ≥ f (x)). Then, if Q(γϵ) > M1, the inequality f (Q(γϵ)) ≥ f (Q(1 - ϵ)
) also holds. In other words, even though a distribution following the mean-γ-median-mode inequality may not be strictly γ-ordered, the inequality defining the γ-orderliness remains valid for most quantile averages. The mean-γ-medianmode inequality can also indicate possible bounds for γ in practice, e.g., for any distributions, when γ → ∞, the γmedian will be greater than the mean and the mode, when γ → 0, the γ-median will be smaller than the mean and the mode, a reasonable γ should maintain the validity of the mean-γ-median-mode inequality.

The definition above of γ-orderliness for a right-skewed distribution implies a monotonic decreasing behavior of the quantile average function with respect to the breakdown point. Therefore, consider the sign of the partial derivative, it can also be expressed as:

∀0 ≤ ϵ ≤ 1 1 + γ , ∂QA ∂ϵ ≤ 0.
The left-skewed case can be obtained by reversing the inequality ∂QA ∂ϵ ≤ 0 to ∂QA ∂ϵ ≥ 0 and employing the second definition of QA, as given in [START_REF] Li | Worst-case range value-at-risk with partial information[END_REF]. For simplicity, the left-skewed case will be omitted in the following discussion. If γ = 1, the γ-ordered distribution is referred to as ordered distribution.

Furthermore, many common right-skewed distributions, such as the Weibull, gamma, lognormal, and Pareto distributions, are partially bounded, indicating a convex behavior of the QA function with respect to ϵ as ϵ approaches 0. By further assuming convexity, the second γ-orderliness can be defined for a right-skewed distribution as follows,

∀0 ≤ ϵ ≤ 1 1 + γ , ∂ 2 QA ∂ϵ 2 ≥ 0 ∧ ∂QA ∂ϵ ≤ 0.
Analogously, the νth γ-orderliness of a right-skewed distribution can be defined as (-1)

ν ∂ ν QA ∂ϵ ν ≥ 0 ∧ . . . ∧ -∂QA ∂ϵ ≥ 0. If γ = 1,
the νth γ-orderliness is referred as to νth orderliness. Let PO denote the set of all distributions that are ordered and PO ν and PγO ν represent the sets of all distributions that are νth ordered and νth γ-ordered, respectively. When the shape parameter of the Weibull distribution, α, is smaller than

D R A F T

cases when it is near-symmetric, as shown in the SI Text. Nevertheless, computing the derivatives of the QA function is often intricate and, at times, challenging. The following theorems establish the relationship between PO, PO ν , and PγO ν , and a wide range of other semi-parametric distributions. They can be used to quickly identify some parametric distributions in PO, PO ν , and PγO ν .

Theorem .2. For any random variable X whose probability distribution function belongs to a location-scale family, the distribution is νth γ-ordered if and only if the family of probability distributions is νth γ-ordered.

Proof. Let Q0 denote the quantile function of the standard distribution without any shifts or scaling. After a locationscale transformation, the quantile function becomes Q(p) = λQ0(p)+µ, where λ is the scale parameter and µ is the location parameter. According to the definition of the νth γ-orderliness, the signs of derivatives of the QA function are invariant after this transformation. As the location-scale transformation is reversible, the proof is complete.

Theorem .2 demonstrates that in the analytical proof of the νth γ-orderliness of a parametric distribution, both the location and scale parameters can be regarded as constants. It is also instrumental in proving other theorems.

Theorem .3. Define a γ-symmetric distribution as one for which the quantile function satisfies

Q(γϵ) = 2Q( γ 1+γ )-Q(1-ϵ) for all 0 ≤ ϵ ≤ 1 1+γ . Any γ-symmetric distribution is νth γ- ordered. Proof. The equality, Q(γϵ) = 2Q( γ 1+γ ) -Q(1 -ϵ), implies that ∂Q(γϵ) ∂ϵ = γQ ′ (γϵ) = ∂(-Q(1-ϵ)) ∂ϵ = Q ′ (1 -ϵ).
From the first definition of QA, the QA function of the γ-symmetric distribution is a horizontal line, since ∂QA ∂ϵ = γQ ′ (γϵ) -Q ′ (1ϵ) = 0. So, the νth order derivative of QA is always zero.

Theorem .4. A symmetric distribution is a special case of the γ-symmetric distribution when

γ = 1, provided that the cdf is monotonic. Proof. A symmetric distribution is a probability distribution such that for all x, f (x) = f (2m -x). Its cdf satisfies F (x) = 1 -F (2m -x). Let x = Q(p), then, F (Q(p)) = p = 1 - F (2m -Q(p)) and F (Q(1 -p)) = 1 -p ⇔ p = 1 -F (Q(1 -p)). Therefore, F (2m -Q(p)) = F (Q(1 -p)). Since the cdf is monotonic, 2m -Q(p) = Q(1 -p) ⇔ Q(p) = 2m -Q(1 -p).
Choosing p = ϵ yields the desired result.

Since the generalized Gaussian distribution is symmetric around the median, it is νth ordered, as a consequence of Theorem .3. Also, the integral of all quantile averages is not equal to the mean, unless γ = 1, as the left and right parts have different weights. The symmetric distribution has a unique role in that its all quantile averages are equal to the mean for a distribution with a finite mean.

Theorem .5. Any right-skewed distribution whose quantile function

Q satisfies Q (ν) (p) ≥ 0 ∧ . . . Q (i) (p) ≥ 0 . . . ∧ Q (2) (p) ≥ 0, i mod 2 = 0, is νth γ-ordered, provided that 0 ≤ γ ≤ 1. Proof. Since (-1) i ∂ i QA ∂ϵ i = 1 2 ((-γ) i Q i (γϵ)+Q i (1-ϵ)) and 1 ≤ i ≤ ν, when i mod 2 = 0, (-1) i ∂ i QA ∂ϵ i ≥ 0 for all γ ≥ 0. When i mod 2 = 1, if further assuming 0 ≤ γ ≤ 1, (-1) i ∂ i QA ∂ϵ i ≥ 0, since Q (i+1) (p) ≥ 0.
This result makes it straightforward to show that the Pareto distribution follows the νth γ-orderliness, provided that 0 ≤ γ ≤ 1, since the quantile function of the Pareto distribution is QP ar (p) = xm(1 -p) -1 α , where xm > 0, α > 0, and so Q (ν) P ar (p) ≥ 0 for all ν ∈ N according to the chain rule. Theorem .6. A right-skewed distribution with a monotonic decreasing pdf is second γ-ordered.

Proof. Given that a monotonic decreasing pdf implies f ′ (x) = F (2) (x) ≤ 0, let x = Q (F (x)), then by differentiating both sides of the equation twice, one can obtain 0 = Q (2)

(F (x)) (F ′ (x)) 2 + Q ′ (F (x)) F (2) (x) ⇒ Q (2) (F (x)) = -Q ′ (F (x))F (2) (x) (F ′ (x)) 2 ≥ 0, since Q ′ (p) ≥ 0.
Theorem .1 already established the γ-orderliness for all γ ≥ 0, which means ∀0 ≤ ϵ ≤ 1 1+γ , ∂QA ∂ϵ ≤ 0. The desired result is then derived from the proof of Theorem .5, since (-1)

2 ∂ 2 QA ∂ϵ 2
≥ 0 for all γ ≥ 0.

Theorem .6 provides valuable insights into the relation between modality and second γ-orderliness. The conventional definition states that a distribution with a monotonic pdf is still considered unimodal. However, within its supported interval, the mode number is zero. Theorem .1 implies that the number of modes and their magnitudes within a distribution are closely related to the likelihood of γ-orderliness being valid. This is because, for a distribution satisfying the necessary and sufficient condition in Theorem .1, it is already implied that the probability density of the left-hand side of the γ-median is always greater than the corresponding probability density of the right-hand side of the γ-median. So although counterexamples can always be constructed for non-monotonic distributions, the general shape of a γ-ordered distribution should have a single dominant mode. It can be easily established that the gamma distribution is second γordered when α ≤ 1, as the pdf of the gamma distribution

is f (x) = λ -α x α-1 e -x λ Γ(α)
, where x ≥ 0, λ > 0, α > 0, and Γ represents the gamma function. This pdf is a product of two monotonic decreasing functions under constraints. For α > 1, analytical analysis becomes challenging. Numerical results can varify that orderliness is valid if α < 140, the second orderliness is valid if α > 81, and the third orderliness is valid if α < 59 (SI Text). It is instructive to consider that when α → ∞, the gamma distribution converges to a Gaussian distribution with mean µ = αλ and variance σ = αλ 2 . The skewness of the gamma distribution, α+2 √ α(α+1)

, is monotonic

with respect to α, since ∂ μ3 (α) ∂α = -3α-2
2(α(α+1)) 3/2 < 0. When α = 59, μ3(α) = 1.025. Theorefore, similar to the Weibull distribution, the parameters which make these distributions fail to be included in PU ∩ PO ∩ PO 2 ∩ PO 3 also correspond to cases when it is near-symmetric.

Theorem .7. Consider a γ-symmetric random variable X.

Let it be transformed using a function ϕ (x) such that ϕ (2) (x) ≥ 0 over the interval supported, the resulting convex transformed

D R A F T distribution is γ-ordered. Moreover, if the quantile function of X satifies Q (2) (p) ≤ 0, the convex transformed distribution is second γ-ordered. Proof. Let ϕQA(ϵ, γ) = 1 2 (ϕ(Q(γϵ)) + ϕ(Q(1 - ϵ))).
Then, for all 0 ≤ ϵ ≤

1 1+γ , ∂ϕQA ∂ϵ = 1 2 (γϕ ′ (Q (γϵ)) Q ′ (γϵ) -ϕ ′ (Q (1 -ϵ)) Q ′ (1 -ϵ)) = 1 2 γQ ′ (γϵ) (ϕ ′ (Q (γϵ)) -ϕ ′ (Q (1 -ϵ))) ≤ 0, since for a γ- symmetric distribution, Q( 1 1+γ )-Q (γϵ) = Q (1 -ϵ)-Q( 1 1+γ ), differentiating both sides, -γQ ′ (γϵ) = -Q ′ (1 -ϵ), where Q ′ (p) ≥ 0, ϕ (2) (x) ≥ 0. If further differentiating the equality, γ 2 Q (2) (γϵ) = -Q (2) (1 -ϵ). Since ∂ (2) ϕQA ∂ϵ (2) = 1 2 γ 2 ϕ 2 (Q (γϵ)) (Q ′ (γϵ)) 2 + ϕ 2 (Q (1 -ϵ)) (Q ′ (1 -ϵ)) 2 + 1 2 γ 2 ϕ ′ (Q (γϵ)) Q 2 (γϵ) + ϕ ′ (Q (1 -ϵ)) Q 2 (1 -ϵ) = 1 2 ϕ (2) (Q (γϵ)) + ϕ (2) (Q (1 -ϵ)) γ 2 Q ′ (γϵ) 2 + 1 2 (ϕ ′ (Q (γϵ)) -ϕ ′ (Q (1 -ϵ))) γ 2 Q (2) (γϵ) . If Q (2) (p) ≤ 0, for all 0 ≤ ϵ ≤ 1 1+γ , ∂ (2)
ϕQA ∂ϵ (2) ≥ 0.

An application of Theorem .7 is that the lognormal distribution is ordered as it is exponentially transformed from the Gaussian distribution. The quantile function of the Gaussian distribution meets the condition Q

(2)

Gau (p) = -2 √ 2πσe 2erfc -1 (2p) 2 erfc -1 (2p) ≤ 0,
where σ is the standard deviation of the Gaussian distribution and erfc denotes the complementary error function. Thus, the lognormal distribution is second ordered. Numerical results suggest that it is also third ordered, although analytically proving this result is challenging.

Theorem .7 also reveals a relation between convex transformation and orderliness, since ϕ is the non-decreasing convex function in van Zwet's trailblazing work Convex transformations of random variables [START_REF] Wr Van Zwet | Convex Transformations of Random Variables: Nebst Stellingen[END_REF] if adding an additional constraint that ϕ ′ (x) ≥ 0. Consider a near-symmetric distribution S, such that the SQA(ϵ) as a function of ϵ fluctuates from 0 to 1 2 . By definition, S is not ordered. Let s be the pdf of S. Applying the transformation ϕ (x) to S decreases s(QS(ϵ)), and the decrease rate, due to the order, is much smaller for s(QS(1 -ϵ)). As a consequence, as ϕ (2) (x) increases, eventually, after a point, for all 0 ≤ ϵ ≤ 1 1+γ , s(QS(ϵ)) becomes greater than s(QS(1 -ϵ)) even if it was not previously. Thus, the SQA(ϵ) function becomes monotonically decreasing, and S becomes ordered. Accordingly, in a family of distributions that differ by a skewness-increasing transformation in van Zwet's sense, violations of orderliness typically occur only when the distribution is near-symmetric.

Pearson proposed using the 3 times standardized meanmedian difference, 3(µ-m) σ , as a measure of skewness in 1895 [START_REF] Pearson | contributions to the mathematical theory of evolution.-ii. skew variation in homogeneous material[END_REF]. [START_REF] Al Bowley | Elements of statistics[END_REF]

if ∀0 ≤ ϵ1 ≤ ϵ2 ≤ 1 2 , SQA ϵ 1 -m ≥ SQA ϵ 2 -m.
Since m is a constant, the monotonic skewness is equivalent to the orderliness. For a nonordered distribution, the signs of SQA ϵ -m with different breakdown points might be different, implying that some skewness measures indicate left-skewed distribution, while others suggest rightskewed distribution. Although it seems reasonable that such a distribution is likely be generally near-symmetric, counterexamples can be constructed. For example, first consider the Weibull distribution, when α > 1 1-ln(2) , it is near-symmetric and nonordered, the non-monotonicity of the SQA function arises when ϵ is close to 1 2 , but if then replacing the third quartile with one from a right-skewed heavy-tailed distribution leads to a right-skewed, heavy-tailed, and nonordered distribution. Therefore, the validity of robust measures of skewness based on the SQA-median difference is closely related to the orderliness of the distribution.

Remarkably, in 2018, Li, Shao, Wang, Yang (2) proved the bias bound of any quantile for arbitrary continuous distributions with finite second moments. Here, let Pµ,σ denotes the set of continuous distributions whose mean is µ and standard deviation is σ. The bias upper bound of the quantile average for P ∈ Pµ=0,σ=1 is given in the following theorem.

Theorem .8. The bias upper bound of the quantile average for any continuous distribution whose mean is zero and standard deviation is one is

sup P ∈P µ=0,σ=1 QA(ϵ, γ) = 1 2 γϵ 1 -γϵ + 1 -ϵ ϵ ,
where 0 ≤ ϵ ≤ 1 1+γ .

Proof. Since sup P ∈P µ=0,σ=1

1 2 (Q(γϵ) + Q(1 -ϵ)) ≤ 1 2 (sup P ∈P µ=0,σ=1 Q(γϵ) + sup P ∈P µ=0,σ=1 Q(1 -ϵ))
, the assertion follows directly from the Lemma 2.6 in (2).

In 2020, Bernard et al. [START_REF] Bernard | Range value-at-risk bounds for unimodal distributions under partial information[END_REF] further refined these bounds for unimodal distributions and derived the bias bound of the symmetric quantile average. Here, the bias upper bound of the quantile average, 0 ≤ γ < 5, for P ∈ PU ∩ Pµ=0,σ=1 is given as sup

P ∈P U ∩P µ=0,σ=1 QA(ϵ, γ) =    1 2 4 9ϵ -1 + 3γϵ 4-3γϵ 0 ≤ ϵ ≤ 1 6 1 2 3(1-ϵ) 4-3(1-ϵ) + 3γϵ 4-3γϵ 1 6 < ϵ ≤ 1 1+γ .
The proof based on the bias bounds of any quantile (3) and the γ ≥ 5 case are given in the SI Text. Subsequent theorems reveal the safeguarding role these bounds play in defining estimators based on νth γ-orderliness.

Theorem .9. sup P ∈P µ=0,σ=1 QA(ϵ, γ) is monotonic decreasing with respect to ϵ over [0, 1 1+γ ], provided that 0 ≤ γ ≤ 1.

Proof. ∂ sup QA(ϵ,γ) ∂ϵ = 1 4 γ γϵ 1-γϵ (γϵ-1) 2 - 1 √ 1 ϵ -1ϵ 2 . When γ = 0, ∂ sup QA(ϵ,γ) ∂ϵ = 1 4 √ γ √ ϵ 1-γϵ (γϵ-1) 2 - 1 √ 1 ϵ -1ϵ 2 = - 1 √ 1 ϵ -1ϵ 2 ≤ 0. When ϵ → 0 + , lim ϵ→0 + 1 4 γ γϵ 1-γϵ (γϵ-1) 2 - 1 √ 1 ϵ -1ϵ 2 = lim ϵ→0 + 1 4 √ γ √ ϵ -1 √ ϵ 3 → -∞. Assuming ϵ > 0, when 0 < γ ≤ 1, to prove ∂ sup QA(ϵ,γ) ∂ϵ ≤ 0, it is equivalent to showing γϵ 1-γϵ (γϵ-1) 2 γ ≥ 1 ϵ -1ϵ 2 . De- fine L(ϵ, γ) = γϵ 1-γϵ (γϵ-1) 2 γ , R(ϵ, γ) = 1 ϵ -1ϵ 2 . Li PNAS | February 11, 2024 | vol. XXX | no. XX | 5 D R A F T L(ϵ,γ) ϵ 2 = γϵ 1-γϵ (γϵ-1) 2 γϵ 2 = 1 γ 1 1 γϵ -1 γ -1 ϵ 2 , R(ϵ,γ) ϵ 2 = 1 ϵ -1. Then, L(ϵ,γ) ϵ 2 ≥ R(ϵ,γ) ϵ 2 ⇔ 1 γ 1 1 γϵ -1 γ -1 ϵ 2 ≥ 1 ϵ -1 ⇔ 1 γ γ -1 ϵ 2 ≥ 1 ϵ -1 1 γϵ -1. Let LmR 1 ϵ = 1 γ 2 γ -1 ϵ 4 -1 ϵ -1 1 γϵ -1 . ∂LmR(1/ϵ) ∂(1/ϵ) = - 4(γ-1 ϵ ) 3 γ 2 - 1 ϵ -1 γ -1 γϵ +1 = -4γ 3 +γ 2 +γ+4 1 ϵ 3 -12γ 1 ϵ 2 +12γ 2 1 ϵ -2γ 1 ϵ γ 2 . Since 0 ≤ γ ≤ 1, 0 ≤ ϵ ≤ 1 1+γ ⇔ 0 ≤ γ ≤ 1 ϵ -1 ⇔ 1 -1 ϵ ≤ -γ ≤ 0 ⇔ 1 ≤ 1 ϵ -γ ≤ 1 ϵ . The numerator of ∂LmR(1/ϵ) ∂(1/ϵ) can be simplified as -4γ 3 + γ 2 + γ + 4 1 ϵ 3 -12γ 1 ϵ 2 + 12γ 2 1 ϵ -2γ 1 ϵ = 4 1 ϵ -γ 3 +γ 2 +γ -2γ 1 ϵ = 4 1 ϵ -γ 3 -γ 2 +γ -2γ 1 ϵ -γ = γ (1 -γ) + 2 1 ϵ -γ 2 1 ϵ -γ 2 -γ . Since 2 1 ϵ -γ 2 ≥ 2, 2 1 ϵ -γ 2 -γ ≥ 2. Also, γ (1 -γ) ≥ 0, 1 ϵ -γ ≥ 0, therefore, γ (1 -γ) + 2 1 ϵ -γ 2 1 ϵ -γ 2 -γ ≥ 0, ∂LmR(1/ϵ) ∂(1/ϵ) ≥ 0. Also, LmR (1 + γ) = 1 γ 2 (γ -1 -γ) 4 - (1 + γ -1) 1 γ (1 + γ) -1 = 1 γ 2 ≥ 0. Therefore, LmR 1 ϵ ≥ 0 for ϵ ∈ (0, 1 1+γ ], provided that 0 < γ ≤ 1.
Consequently, the simplified inequality

1 γ γ -1 ϵ 2 ≥ 1 ϵ -1 1 γϵ -1 is valid. ∂ sup QA(ϵ,γ) ∂ϵ is non-positive throughout the interval 0 ≤ ϵ ≤ 1 1+γ , given that 0 ≤ γ ≤ 1, the proof is complete. Theorem .10. sup P ∈P U ∩P µ=0,σ=1 QA(ϵ, γ) is a nonincreasing function with respect to ϵ on the interval [0, 1 1+γ ], provided that 0 ≤ γ ≤ 1. Proof. When 0 ≤ ϵ ≤ 1 6 , ∂ sup QA ∂ϵ = γ ϵγ 12-9ϵγ (4-3ϵγ) 2 - 1 3 √ 4 ϵ -9ϵ 2 = √ γ √ ϵ 12-9ϵγ (4-3ϵγ) 2 - 1 3 √ 4 ϵ -9ϵ 2 . If γ = 0 and ϵ → 0 + , ∂ sup QA ∂ϵ = - 1 3 √ 4 ϵ -9ϵ 2 < 0. If ϵ → 0 + , lim ϵ→0 + γ (4-3γϵ) 2 ϵγ 12-9γϵ - 1 3 √ 4 ϵ -9ϵ 2 = lim ϵ→0 + √ 3γ √ 4 3 ϵ -1 6 √ ϵ 3 → -∞, for all 0 ≤ γ ≤ 1, so, ∂ sup QA ∂ϵ < 0. When 0 < ϵ ≤ 1 6 and 0 < γ ≤ 1, to prove ∂ sup QA ∂ϵ ≤ 0, it is equivalent to showing ϵγ 12-9ϵγ (4-3ϵγ) 2 γ ≥ 3 4 ϵ -9ϵ 2 . Define L(ϵ, γ) = ϵγ 12-9ϵγ (4-3ϵγ) 2 γ , R(ϵ, γ) = 3 4 ϵ -9ϵ 2 . L(ϵ,γ) ϵ 2 = ϵγ 12-9ϵγ (4-3ϵγ) 2 γϵ 2 = 1 γ 4 ϵ -3γ 2 1 12 ϵγ -9 , R(ϵ,γ) ϵ 2 = 3 4 ϵ -9.
Then, the above inequality is

equivalent to L(ϵ,γ) ϵ 2 ≥ R(ϵ,γ) ϵ 2 ⇔ 1 γ 1 12 ϵγ -9 4 ϵ -3γ 2 ≥ 3 4 ϵ -9 ⇔ 1 γ 4 ϵ -3γ 2 ≥ 3 12 ϵγ -9 4 ϵ -9 ⇔ 1 γ 2 4 ϵ -3γ 4 ≥ 9 12 ϵγ -9 4 ϵ -9 . Let LmR 1 ϵ = 1 γ 2 4 ϵ -3γ 4 -9 12 ϵγ -9 4 ϵ -9 . ∂LmR(1/ϵ) ∂(1/ϵ) = 16( 4 ϵ -3γ) 3 γ 2 - 36 12 ϵγ -9 - 108(4 4 ϵ -9) γ = 4(4( 4 ϵ -3γ) 3 -27γ( 4 ϵ -3γ)+27(9-4 ϵ )γ) γ 2 = 4(256 1 ϵ 3 -576 1 ϵ 2 γ+432 1 ϵ γ 2 -216 1 ϵ γ-108γ 3 +81γ 2 +243γ) γ 2 .
Since

256 1 ϵ 3 -576 1 ϵ 2 γ + 432 1 ϵ γ 2 -216 1 ϵ γ -108γ 3 + 81γ 2 + 243γ ≥ 1536 1 ϵ 2 -576 1 ϵ 2 + 432 1 ϵ γ 2 -216 1 ϵ γ -108γ 3 + 81γ 2 + 243γ ≥ 924 1 ϵ 2 + 36 1 ϵ 2 -216 1 ϵ + 432 1 ϵ γ 2 -108γ 3 + 81γ 2 + 243γ ≥ 924 1 ϵ 2 + 36 1 ϵ 2 -216 1 ϵ + 513γ 2 -108γ 3 + 243γ > 0, ∂LmR(1/ϵ) ∂(1/ϵ) > 0. Also, LmR (6) = 81(γ-8)((γ-8) 3 +15γ) γ 2 > 0 ⇐⇒ γ 4 -32γ 3 + 399γ 2 -2168γ + 4096 > 0. If 0 < γ ≤ 1,
then 32γ 3 < 256. Also, γ 4 > 0. So, it suffices to prove that 399γ 2 -2168γ + 4096 > 256. Applying the quadratic formula demonstrates the validity of LmR (6) > 0, if 0 < γ ≤ 1. Hence, LmR 1 ϵ ≥ 0 for ϵ ∈ (0, 1 6 ], if 0 < γ ≤ 1. The first part is finished. When

1 6 < ϵ ≤ 1 1+γ , ∂ sup QA ∂ϵ = √ 3 γ √ γϵ(4-3γϵ) 3 2 - 1 √ 1-ϵ(3ϵ+1) 3 2
. If γ = 0,

γ √ γϵ(4-3γϵ) 3 2 = √ γ √ ϵ(4-3γϵ) 3 2 = 0, so ∂ sup QA ∂ϵ = √ 3 - 1 √ 1-ϵ(3ϵ+1) 3 2 < 0, for all 1 6 < ϵ ≤ 1 1+γ . If γ > 0, to determine whether ∂ sup QA ∂ϵ ≤ 0, when 1 6 < ϵ ≤ 1 1+γ , since √ 1 -ϵ (3ϵ + 1) 3 2 > 0 and √ γϵ (4 -3γϵ) 3 2 > 0, showing √ γϵ(4-3γϵ) 3 2 γ ≥ √ 1 -ϵ (3ϵ + 1) 3 2 ⇔ γϵ(4-3γϵ) 3 γ 2 ≥ (1 -ϵ) (3ϵ + 1) 3 ⇔ -27γ 2 ϵ 4 + 108γϵ 3 + 64ϵ γ + 27ϵ 4 -162ϵ 2 -8ϵ -1 ≥ 0 is sufficient. When 0 < γ ≤ 1, the inequality can be further simplified to 108γϵ 3 + 64ϵ γ -162ϵ 2 -8ϵ -1 ≥ 0. Since ϵ ≤ 1 1+γ , γ ≤ 1 ϵ -1. Also, as 0 < γ ≤ 1 is assumed, the range of γ can be expressed as 0 < γ ≤ min(1, 1 ϵ -1). When 1 6 < ϵ ≤ 1 2 , 1 < 1 ϵ -1, so in this case, 0 < γ ≤ 1. When 1 2 ≤ ϵ < 1, so in this case, 0 < γ ≤ 1 ϵ -1. Let h(γ) = 108γϵ 3 + 64ϵ γ , ∂h(γ) ∂γ = 108ϵ 3 -64ϵ γ 2 . When γ ≤ 64ϵ 18ϵ 3 , ∂h(γ) ∂γ ≥ 0, when γ ≥ 64ϵ 18ϵ 3 , ∂h(γ)
∂γ ≤ 0, therefore, the minimum of h(γ) must be when γ is equal to the boundary point of the domain. When 1 6 < ϵ ≤ 1 2 , 0 < γ ≤ 1, since h(0) → ∞, h(1) = 108ϵ 3 +64ϵ, the minimum occurs at the boundary point QA with respect to ϵ implies that the extent of any violations of the γ-orderliness, if 0 ≤ γ ≤ 1, is bounded for any distribution with a finite second moment, e.g., for a right-skewed distribution in

γ = 1, 108γϵ 3 + 64ϵ γ -162ϵ 2 -8ϵ-1 > 108ϵ 3 +56ϵ-162ϵ 2 -1. Let g(ϵ) = 108ϵ 3 + 56ϵ -162ϵ 2 -1. g ′ (ϵ) = 324ϵ 2 -324ϵ + 56, when ϵ ≤ 2 9 , g ′ (ϵ) ≥ 0, when 2 9 ≤ ϵ ≤ 1 2 , g ′ (ϵ) ≤ 0, since g( 1 6 ) = 13 3 , g( 1 2 ) = 0, so g(ϵ) ≥ 0, 108γϵ 3 + 64ϵ γ -162ϵ 2 -8ϵ -1 ≥ 0. When 1 2 ≤ ϵ < 1, 0 < γ ≤ 1 ϵ -1. Since h( 1 ϵ -1) = 108( 1 ϵ -1)ϵ 3 + 64ϵ 1 ϵ -1 , 108γϵ 3 + 64ϵ γ -162ϵ 2 -8ϵ -1 > 108 1 ϵ -1 ϵ 3 + 64ϵ 1 ϵ -1 -162ϵ 2 -8ϵ -1 = -108ϵ 4 +54ϵ 3 -18ϵ 2 +7ϵ+1 ϵ-1 . Let nu(ϵ) = -108ϵ 4 + 54ϵ 3 -18ϵ 2 + 7ϵ + 1, then nu ′ (ϵ) = -432ϵ 3 + 162ϵ 2 -36ϵ + 7, nu ′′ (ϵ) = -1296ϵ 2 + 324ϵ -36 < 0. Since nu ′ (ϵ = 1 2 ) = -49 2 < 0, nu ′ (ϵ) < 0. Also, nu(ϵ = 1 2 ) = 0, so nu(ϵ) ≥ 0, 108γϵ 3 + 64ϵ γ -162ϵ 2 -8ϵ -1 ≥ 0 is also valid. As a result,
P 2 Υ , if 0 ≤ ϵ1 ≤ ϵ2 ≤ ϵ3 ≤ 1 1+γ , QA ϵ 2 ,γ ≥ QA ϵ 3 ,γ ≥ QA ϵ 1 ,
γ , then QA ϵ 2 ,γ will not be too far away from QA ϵ 1 ,γ , since sup

P ∈P 2 Υ QA ϵ 1 ,γ > sup P ∈P 2 Υ QA ϵ 2 ,γ > sup P ∈P 2 Υ QA ϵ 3 ,γ .
Moreover, a stricter bound can be established for unimodal distributions according to Bernard et al. 's result [START_REF] Bernard | Range value-at-risk bounds for unimodal distributions under partial information[END_REF]. The violation of νth γ-orderliness, when ν ≥ 2, is also bounded, since the QA function is bounded, the νth γ-orderliness corresponds to the higher-order derivatives of the QA function with respect to ϵ.

The Impact of γ-Orderliness on Weighted Inequalities

Analogous to the γ-orderliness, the γ-trimming inequality for a right-skewed distribution is defined as ∀0 ≤ ϵ1 ≤ ϵ2 ≤ 1 1+γ , TMϵ 1 ,γ ≥ TMϵ 2 ,γ . γ-orderliness is a sufficient condition for the γ-trimming inequality, as proven in the SI Text. The next theorem shows a relation between the ϵ,γ-quantile average and the ϵ,γ-trimmed mean under the γ-trimming inequality, suggesting the γ-orderliness is not a necessary condition for the γ-trimming inequality.

Theorem .11. For a distribution that is right-skewed and follows the γ-trimming inequality, it is asymptotically true that the quantile average is always greater or equal to the corresponding trimmed mean with the same ϵ and γ, for all 0 ≤ ϵ ≤ 1 1+γ . Proof. According to the definition of the γ-trimming inequality:

∀0 ≤ ϵ ≤ 1 1+γ , 1 1-ϵ-γϵ+2δ 1-ϵ+δ γϵ-δ Q (u) du ≥ 1 1-ϵ-γϵ 1-ϵ γϵ Q (u) du,
where δ is an infinitesimal positive quantity.

Subsequently, rewriting the inequality gives

1-ϵ+δ γϵ-δ Q (u) du -1-ϵ-γϵ+2δ 1-ϵ-γϵ 1-ϵ γϵ Q (u) du ≥ 0 ⇔ 1-ϵ+δ 1-ϵ Q (u) du + γϵ γϵ-δ Q (u) du - 2δ 1-ϵ-γϵ 1-ϵ γϵ Q (u) du ≥ 0. Since δ → 0 + , 1 2δ 1-ϵ+δ 1-ϵ Q (u) du + γϵ γϵ-δ Q (u) du = Q(γϵ)+Q(1-ϵ) 2 ≥ 1 1-ϵ-γϵ 1-ϵ γϵ Q (u) du, the proof is com- plete.
An analogous result about the relation between the ϵ,γtrimmed mean and the ϵ,γ-Winsorized mean can be obtained in the following theorem.

Theorem .12. For a right-skewed distribution following the γ-trimming inequality, asymptotically, the Winsorized mean is always greater or equal to the corresponding trimmed mean with the same ϵ and γ, for all 0 ≤ ϵ ≤ 1 1+γ , provided that 0 ≤ γ ≤ 1. If assuming γ-orderliness, the inequality is valid for any non-negative γ.

Proof. According to Theorem .11, Replacing the TM in the γ-trimming inequality with WA forms the definition of the γ-weighted inequality. The γorderliness also implies the γ-Winsorization inequality when 0 ≤ γ ≤ 1, as proven in the SI Text. The same rationale as presented in Theorem .2, for a location-scale distribution characterized by a location parameter µ and a scale parameter λ, asymptotically, any WA(ϵ, γ) can be expressed as λWA0(ϵ, γ) + µ, where WA0(ϵ, γ) is an function of Q0(p) according to the definition of the weighted average. Adhering to the rationale present in Theorem .2, for any probability distribution within a location-scale family, a necessary and sufficient condition for whether it follows the γ-weighted inequality is whether the family of probability distributions also adheres to the γ-weighted inequality.

Q(γϵ)+Q(1-ϵ) 2 ≥ 1 1-ϵ-γϵ 1-ϵ γϵ Q (u) du ⇔ γϵ (Q (γϵ) + Q (1 -ϵ)) ≥ ( 2γϵ 1-ϵ-γϵ ) 1-ϵ γϵ Q (u) du. Then, if 0 ≤ γ ≤ 1, 1 - 1 1-ϵ-γϵ 1-ϵ γϵ Q (u) du + γϵ (Q (γϵ) + Q (1 -ϵ)) ≥ 0 ⇒ 1-ϵ γϵ Q (u) du + γϵQ (γϵ) + ϵQ (1 -ϵ) ≥ 1-ϵ γϵ Q (u) du + γϵ (Q (γϵ) + Q (1 -ϵ)) ≥ 1 1-ϵ-γϵ 1-ϵ γϵ Q (u)
To construct weighted averages based on the νth γorderliness and satisfying the corresponding weighted inequality, when 0

≤ γ ≤ 1, let Bi = (i+1)ϵ iϵ QA (u, γ) du, ka = kϵ + c. From the γ-orderliness for a right-skewed dis- tribution, it follows that, -∂QA ∂ϵ ≥ 0 ⇔ ∀0 ≤ a ≤ 2a ≤ 1 1+γ , -(QA(2a,γ)-QA(a,γ)) a ≥ 0 ⇒ Bi -Bi+1 ≥ 0, if 0 ≤ γ ≤ 1.
Suppose that Bi = B0. Then, the ϵ,γ-block Winsorized mean, is defined as

BWMϵ,γ,n := 1 n   (1-ϵ)n i=nγϵ+1 Xi + 2nγϵ+1 i=nγϵ+1 Xi + (1-ϵ)n i=(1-2ϵ)n Xi   ,
which is double weighting the leftest and rightest blocks having sizes of γϵn and ϵn, respectively. As a consequence of Bi -Bi+1 ≥ 0, the γ-block Winsorization inequality is valid, provided that 0 ≤ γ ≤ 1. The block Winsorized mean uses two blocks to replace the trimmed parts, not two single quantiles. The subsequent theorem provides an explanation for this difference.

Theorem .13. Asymptotically, for a right-skewed distribution following the γ-orderliness, the Winsorized mean is always greater than or equal to the corresponding block Winsorized mean with the same ϵ and γ, for all 0 ≤ ϵ ≤ 1 1+γ , provided that 0 ≤ γ ≤ 1.

Proof. From the definitions of BWM and WM, the statement necessitates

1-ϵ γϵ Q (u) du + γϵQ (γϵ) + ϵQ (1 -ϵ) ≥ 1-ϵ γϵ Q (u) du + 2γϵ γϵ Q (u) du + 1-ϵ 1-2ϵ Q (u) du ⇔ γϵQ (γϵ) + ϵQ (1 -ϵ) ≥ 2γϵ γϵ Q (u) du+ 1-ϵ 1-2ϵ Q (u) du. Define WMl(x) = Q (γϵ) and BWMl(x) = Q (x).
In both functions, the interval for x is specified as [γϵ, 2γϵ].

Then, define WMu(y) = Q (1 -ϵ) and BWMu(y) = Q (y). In both functions, the interval for y is specified as

[1 -2ϵ, 1 -ϵ]. The function y : [γϵ, 2γϵ] → [1 -2ϵ, 1 -ϵ] defined by y(x) = 1 -x γ is a bijection. WMl(x) + WMu(y(x)) = Q (γϵ) + Q (1 -ϵ) ≥ BWMl(x) + BWMu(y(x)) = Q (x) + Q 1 -x γ
is valid for all x ∈ [γϵ, 2γϵ], according to the definition of γ-orderliness. Integration of the left side yields,

2γϵ γϵ (WMl (u) + WMu (y (u))) du = 2γϵ γϵ Q (γϵ) du + y(2γϵ) y(γϵ) Q (1 -ϵ) du = 2γϵ γϵ Q (γϵ) du + 1-ϵ 1-2ϵ Q (1 -ϵ) du = γϵQ (γϵ) + ϵQ (1 -ϵ), while integration of the right side yields 2γϵ γϵ (BWMl (x) + BWMu (y (x))) dx = 2γϵ γϵ Q (u) du + 2γϵ γϵ Q 1 -x γ dx = 2γϵ γϵ Q (u) du + 1-ϵ 1-2ϵ Q (u) du,

D R A F T

upper limits and lower limits of the integrations are different for each term, the condition 0 ≤ γ ≤ 1 is necessary for the desired inequality to be valid.

From the second γ-orderliness for a right-skewed distribution,

∂ 2 QA ∂ 2 ϵ ≥ 0 ⇒ ∀0 ≤ a ≤ 2a ≤ 3a ≤ 1 1+γ , 1 a (QA(3a,γ)-QA(2a,γ)) a -(QA(2a,γ)-QA(a,γ)) a ≥ 0 ⇒ if 0 ≤ γ ≤ 1, Bi -2Bi+1 + Bi+2 ≥ 0.
SMϵ can thus be interpreted as assuming γ = 1 and replacing the two blocks, Bi + Bi+2 with one block 2Bi+1. From the νth γ-orderliness for a rightskewed distribution, the recurrence relation of the derivatives naturally produces the alternating binomial coefficients,

(-1) ν ∂ ν QA ∂ϵ ν ≥ 0 ⇒ ∀0 ≤ a ≤ . . . ≤ (ν + 1)a ≤ 1 1 + γ , (-1) ν a   QA(νa+a,γ) . . . a - . . . QA(2a,γ) a a - QA(νa,γ) . . . a - . . . QA(a,γ) a a   ≥ 0 ⇔ (-1) ν a ν ν j=0 (-1) j ν j QA ((ν -j + 1) a, γ) ≥ 0 ⇒ if 0 ≤ γ ≤ 1, ν j=0 (-1) j ν j Bi+j ≥ 0.
Based on the νth orderliness, the ϵ,γ-binomial mean is introduced as

BMν,ϵ,γ,n := 1 n   1 2 ϵ -1 (ν+1) -1 i=1 ν j=0 1 -(-1) j ν j Bi j   ,
where Bi j = nϵ(j+(i-1)(ν+1)+1) l=nγϵ(j+(i-1)(ν+1))+1 (X l + X n-l+1 ). If ν is not indicated, it defaults to ν = 3. Since the alternating sum of binomial coefficients equals zero, when ν ≪ ϵ -1 and ϵ → 0, BM → µ. The solutions for the continuity of the breakdown point is the same as that in SM and not repeated here. The equalities BMν=1,ϵ = BWMϵ and BMν=2,ϵ = SM ϵ,b=3 hold, when γ = 1 and their respective ϵs are identical. Interestingly, the biases of the SM ϵ= 1 9 ,b=3 and the WM ϵ= 1 9 are nearly indistinguishable in common asymmetric unimodal distributions such as Weibull, gamma, lognormal, and Pareto (SI Dataset S1). This indicates that their robustness to departures from the symmetry assumption is practically similar under unimodality, even though they are based on different orders of orderliness. If single quantiles are used, based on the second γ-orderliness, the stratified quantile mean can be defined as SQM ϵ,γ,n := 4ϵ

1 4ϵ i=1 1 2 ( Qn ((2i -1)γϵ) + Qn (1 -(2i -1)ϵ)), SQM ϵ= 1 4
is the Tukey's midhinge [START_REF] Jw Tukey | Exploratory data analysis[END_REF]. In fact, SQM is a subcase of SM when γ = 1 and b → ∞, so the solution for the continuity of the breakdown point, 1 ϵ mod 4 ̸ = 0, is identical. However, since the definition is based on the empirical quantile function, no decimal issues related to order statistics will arise. The next theorem explains another advantage.

Theorem .14. For a right-skewed second γ-ordered distribution, asymptotically, SQM ϵ,γ is always greater or equal to the corresponding BMν=2,ϵ,γ with the same ϵ and γ, for all

0 ≤ ϵ ≤ 1 1+γ , if 0 ≤ γ ≤ 1.
Proof. For simplicity, suppose the order statistics of the sample are distributed into ϵ -1 ∈ N blocks in the computation of both SQM ϵ,γ and BMν=2,ϵ,γ. The computation of BMν=2,ϵ,γ alternates between weighting and non-weighting, let '0' denote the block assigned with a weight of zero and '1' denote the block assigned with a weighted of one, the sequence indicating the weighted or non-weighted status of each block is: 0, 1, 0, 0, 1, 0, . . . . Let this sequence be denoted by a BM ν=2,ϵ,γ (j), its formula is a BM ν=2,ϵ,γ (j) = j mod 3

2

. Similarly, the computation of SQM ϵ,γ can be seen as positioning quantiles (p) at the beginning of the blocks if 0 < p < 1 1+γ , and at the end of the blocks if p > 1 1+γ . The sequence of denoting whether each block's quantile is weighted or not weighted is: 0, 1, 0, 1, 0, 1, . . . . Let the sequence be denoted by a SQM ϵ,γ (j), the formula of the sequence is a SQM ϵ,γ (j) = j mod 2. If pairing all blocks in BMν=2,ϵ,γ and all quantiles in SQM ϵ,γ , there are two possible pairings of a BM ν=2 (j) and a SQM ϵ,γ (j). One pairing occurs when a BM ν=2,ϵ,γ (j) = a SQM ϵ,γ (j) = 1, while the other involves the sequence 0, 1, 0 from a BM ν=2,ϵ,γ (j) paired with 1, 0, 1 from a SQM ϵ,γ (j). By leveraging the same principle as Theorem .13 and the second γ-orderliness (replacing the two quantile averages with one quantile average between them), the desired result follows.

The biases of SQM ϵ= 1 8 , which is based on the second orderliness with a quantile approach, are notably similar to those of BM ν=3,ϵ= 1 8 , which is based on the third orderliness with a block approach, in common asymmetric unimodal distributions (Figure 1).

Hodges-Lehmann Inequality and γ-U -Orderliness

The Hodges-Lehmann estimator stands out as a unique robust location estimator due to its definition being substantially dissimilar from conventional L-estimators, R-estimators, and M -estimators. In their landmark paper, Estimates of location based on rank tests, Hodges and Lehmann [START_REF] Hodges | Estimates of location based on rank tests[END_REF] proposed two methods for computing the H-L estimator: the Wilcoxon score R-estimator and the median of pairwise means. The Wilcoxon score R-estimator is a location estimator based on signedrank test, or R-estimator, [START_REF] Hodges | Estimates of location based on rank tests[END_REF] and was later independently discovered by Sen (1963) [START_REF] Pk Sen | On the estimation of relative potency in dilution (-direct) assays by distribution-free methods[END_REF]. However, the median of pairwise means is a generalized L-statistic and a trimmed U -statistic, as classified by Serfling in his novel conceptualized study in 1984 [START_REF] Serfling | Generalized l-, m-, and r-statistics[END_REF]. Serfling further advanced the understanding by generalizing the H-L kernel as

hl k (x1, . . . , x k ) = 1 k k i=1 xi,
where k ∈ N (45). Here, the weighted H-L kernel is defined

as whl k (x1, . . . , x k ) = k i=1 x i w i k i=1 w i
, where wis are the weights applied to each element.

By using the weighted H-L kernel and the L-estimator, it is now clear that the Hodges-Lehmann estimator is an LLstatistic, the definition of which is provided as follows:

LL k,ϵ,γ,n := Lϵ 0 ,γ,n sort (whl k (XN 1 , ••• , XN k )) ( n k ) N =1
, where Lϵ 0 ,γ,n (Y ) represents the ϵ0,γ-L-estimator that uses the sorted sequence, sort (whl

k (XN 1 , ••• , XN k )) ( n k ) N =1
, as input. The upper asymptotic breakdown point of LL k,ϵ,γ is ϵ = 1 -(1 -ϵ0) 1 k , as proven in REDS III [START_REF] Li | Robust estimations from distribution structures: Iii. invariant moments[END_REF]. There are two ways to adjust the breakdown point: either by setting k as a

| Li

D R A F T

constant and adjusting ϵ0, or by setting ϵ0 as a constant and adjusting k. In the above definition, k is discrete, but the bootstrap method can be applied to ensure the continuity of k, also making the breakdown point continuous. Specifically, if k ∈ R, let the bootstrap size be denoted by b, then first sampling the original sample (1 -k + ⌊k⌋)b times with each sample size of ⌊k⌋, and then subsequently sampling (1 -⌈k⌉ + k)b times with each sample size of ⌈k⌉, (1

-k + ⌊k⌋)b ∈ N, (1 -⌈k⌉ + k)b ∈ N.
The corresponding kernels are computed separately, and the pooled sorted sequence is used as the input for the L-estimator.

Let S k represent the sorted sequence. Indeed, for any finite sample, X, when k = n, S k becomes a single point, whl k=n (X1, . . . , Xn). When wi = 1, the minimum of

S k is 1 k k i=1
Xi, due to the property of order statistics. The maximum of S k is 1 k k i=1 Xn-i+1. The monotonicity of the order statistics implies the monotonicity of the extrema with respect to k, i.e., the support of S k shrinks monotonically. For unequal wis, the shrinkage of the support of S k might not be strictly monotonic, but the general trend remains, since all LL-statistics converge to the same point, as k → n. Therefore,

if n i=1 X i w i n i=1 w i
approaches the population mean when n → ∞, all LL-statistics based on such consistent kernel function approach the population mean as k → ∞. For example, if

whl k = BM ν,ϵ k ,n=k , ν ≪ ϵ -1
k , ϵ k → 0, such kernel function is consistent. These cases are termed the LL-mean (LLM k,ϵ,γ,n ). By substituting the WAϵ 0 ,γ,n for the Lϵ 0 ,γ,n in LL-statistic, the resulting statistic is referred to as the weighted L-statistic (WL k,ϵ,γ,n ). The case having a consistent kernel function is termed as the weighted L-mean (WLM k,ϵ,γ,n ). The wi = 1 case of WLM k,ϵ,γ,n is termed the weighted Hodges-Lehmann mean (WHLM k,ϵ,γ,n ). The WHLM k=1,ϵ,γ,n is the weighted average. If k ≥ 2 and the WA in WHLM is set as TMϵ 0 , it is called the trimmed H-L mean (Figure 1, k = 2, ϵ0 = 15 64 ). The THLM k=2,ϵ,γ=1,n appears similar to the Wilcoxon's onesample statistic investigated by Saleh in 1976 [START_REF] Saleh | Hodges-lehmann estimate of the location parameter in censored samples[END_REF], which involves first censoring the sample, and then computing the mean of the number of events that the pairwise mean is greater than zero. The THLM

k=2,ϵ=1-(1-1 2 ) 1 2 ,γ=1,n
is the Hodges-Lehmann estimator, or more generally, a special case of the median Hodges-Lehmann mean (mHLM k,n ). mHLM k,n is asymptotically equivalent to the MoM k,b= n k as discussed previously, Therefore, it is possible to define a series of location estimators, analogous to the WHLM, based on MoM. For example, the γ-median of means, γmoM k,b= n k ,n , is defined by replacing the median in MoM k,b= n k ,n with the γ-median.

The hl k kernel distribution, denoted as F hl k , can be defined as the probability distribution of the sorted sequence,

sort (hl k (XN 1 , ••• , XN k )) ( n k ) N =1
. For any real value y, the cdf of the hl k kernel distribution is given by: F h k (y) = P(Yi ≤ y), where Yi represents an individual element from the sorted sequence. The overall hl k kernel distributions possess a twodimensional structure, encompassing n kernel distributions with varying k values, from 1 to n, where one dimension is inherent to each individual kernel distribution, while the other is formed by the alignment of the same percentiles across all kernel distributions. As k increases, all percentiles converge to X, leading to the concept of γ-U -orderliness:

(∀k2 ≥ k1 ≥ 1, γmHLM k 2 ,ϵ=1- γ 1+γ 1 k 2 ,γ ≥ γmHLM k 1 ,ϵ=1- γ 1+γ 1 k 1 ,γ )∨ (∀k2 ≥ k1 ≥ 1, γmHLM k 2 ,ϵ=1- γ 1+γ 1 k 2 ,γ ≤ γmHLM k 1 ,ϵ=1- γ 1+γ 1 k 1 ,γ ),
where γmHLM k sets the WA in WHLM as γ-median, with γ being constant. The direction of the inequality depends on the relative magnitudes of γmHLM k=1,ϵ,γ = γm and γmHLM k=∞,ϵ,γ = µ. The Hodges-Lehmann inequality can be defined as a special case of the γ-U -orderliness when γ = 1. When γ ∈ {0, ∞}, the γ-U -orderliness is valid for any distribution as previously shown. If γ / ∈ {0, ∞}, analytically proving the validity of the γ-U -orderliness for a parametric distribution is pretty challenging. As an example, the hl2 kernel distribution has a probability density function f hl 2 (x) = 2x 0 2f (t) f (2x -t) dt (a result after the transformation of variables); the support of the original distribution is assumed to be [0, ∞) for simplicity. The expected value of the H-L estimator is the positive solution of

H-L 0 (f hl 2 (s)) ds = 1 2 . For the exponential distribution, f hl 2 ,exp (x) = 4λ -2 xe -2λ -1 x , λ is a scale parameter, E[H-L] = -W -1( -1 2e )-1 2 λ ≈ 0.839λ,
where W-1 is a branch of the Lambert W function which cannot be expressed in terms of elementary functions. However, the violation of the γ-U -orderliness is bounded under certain assumptions, as shown below.

Theorem .15. For any distribution with a finite second central moment, σ 2 , the following concentration bound can be established for the γ-median of means,

P γmoM k,b= n k ,n -µ > tσ √ k ≤ e -2n k 1 1+γ -1 k+t 2 2 .
Proof. Denote the mean of each block as µi,

,n -µ > tσ √ k ⊂ b i=1 1 µ i -µ > tσ √ k ≥ b 1 -γ 1+γ
, where 1A is the indicator of event A. Assuming a finite second central moment, σ 2 , it follows from one-sided Chebeshev's inequality that

E 1 µ i -µ > tσ √ k = P ( µi -µ) > tσ √ k ≤ σ 2 kσ 2 +t 2 σ 2 . Given that 1 µ i -µ > tσ √ k ∈ [0, 1] are independent
and identically distributed random variables, according to the aforementioned inclusion relation, the onesided Chebeshev's inequality and the one-sided Hoeffding's inequality,

P γmoM k,b= n k ,n -µ > tσ √ k ≤ P b i=1 1 µ i -µ > tσ √ k ≥ b 1 -γ 1+γ = P 1 b b i=1 1 µ i -µ > tσ √ k -E 1 µ i -µ > tσ √ k ≥ 1 -γ 1+γ -E 1 µ i -µ > tσ √ k ≤ e -2b 1- γ 1+γ -E 1 µ i -µ > tσ √ k 2 ≤ e -2b 1- γ 1+γ - σ 2 kσ 2 +t 2 σ 2 2 = e -2b 1 1+γ -1 k+t 2 2 . Li PNAS | February 11, 2024 | vol. XXX | no. XX | 9 D R A F T Theorem .16. Let B(k, γ, t, n) = e -2n k 1 1+γ -1 k+t 2 2 . If n ∈ N, γ ≥ 0, 0 ≤ t 2 < γ + 1, and γ -t 2 + 1 ≤ k ≤ 1 2 9γ 2 + 18γ -8γt 2 -8t 2 + 9+ 1 2 3γ -2t 2 + 3 , B is mono- tonic decreasing with respect to k. Proof. Since ∂B ∂k = 2n 1 γ+1 -1 k+t 2 2 k 2 - 4n 1 γ+1 -1 k+t 2 k(k+t 2 ) 2 e - 2n 1 γ+1 -1 k+t 2 2 k and n ∈ N, ∂B ∂k ≤ 0 ⇔ 2n 1 γ+1 -1 k+t 2 2 k 2 - 4n 1 γ+1 -1 k+t 2 k(k+t 2 ) 2 ≤ 0 ⇔ 2n(-γ+k+t 2 -1)(k 2 -3(γ+1)k+2kt 2 +t 2 (-γ+t 2 -1)) (γ+1) 2 k 2 (k+t 2 ) 3 ≤ 0 ⇔ -γ + k + t 2 -1 k 2 -3(γ + 1)k + 2kt 2 + t 2 -γ + t 2 -1 ≤ 0.
When the factors are expanded, it yields a cubic inequality in terms of k:

k 3 + k 2 3t 2 -4(γ + 1) + 3k γ -t 2 + 1 2 + t 2 γ -t 2 + 1 2 ≤ 0. Assuming 0 ≤ t 2 < γ + 1 and γ ≥ 0,
using the factored form and subsequently applying the quadratic formula, the inequality is valid if

γ -t 2 + 1 ≤ k ≤ 1 2 9γ 2 + 18γ -8γt 2 -8t 2 + 9 + 1 2 3γ -2t 2 + 3 .
Let X be a random variable and Ȳ = 1 k (Y1 + • • • + Y k ) be the average of k independent, identically distributed copies of X. Applying the variance operation gives:

Var( Ȳ ) = Var 1 k (Y1 + • • • + Y k ) = 1 k 2 (Var(Y1) + • • • + Var(Y k )) = 1 k 2 (kσ 2 ) = σ 2
k , since the variance operation is a linear operator for independent variables, and the variance of a scaled random variable is the square of the scale times the variance of the variable, i.e., Var(cX

) = E[(cX -E[cX]) 2 ] = E[(cX-cE[X]) 2 ] = E[c 2 (X-E[X]) 2 ] = c 2 E[((X)-E[X]) 2 ] = c 2
Var(X). Thus, the standard deviation of the hl k kernel distribution, asymptotically, is σ √ k . By utilizing the asymptotic bias bound of any quantile for any continuous distribution with a finite second central moment, σ 2 (2), a conservative asymptotic bias bound of γmoM k,b= n k can be established as γmoM

k,b= n k -µ ≤ γ 1+γ 1- γ 1+γ σ hl k = γ k σ. That implies in Theorem .15, t < √ γ, so when γ = 1, the upper bound of k, subject to the monotonic decreasing constraint, is 2 + √ 5 < 1 2 √ 9 + 18 -8t 2 -8t 2 + 9 + 1 2 3 -2t 2 + 3 ≤ 6, the lower bound is 1 < 2 -t 2 ≤ 2.
These analyses elucidate a surprising result: although the conservative asymptotic bound of MoM k,b= n k is monotonic with respect to k, its concentration bound is optimal when k ∈ (2 + √ 5, 6]. Then consider the structure within each individual hl k kernel distribution. The sorted sequence S k , when k = n -1, has n elements and the corresponding hl k kernel distribution can be seen as a location-scale transformation of the original distribution, so the corresponding hl k kernel distribution is νth γ-ordered if and only if the original distribution is νth γ-ordered according to Theorem .2. Analytically proving other cases is challenging. For example,

f ′ hl 2 (x) = 4f (2x) f (0)+ 2x 0 4f (t) f ′ (2x -t) dt, the strict neg- ative of f ′ hl 2 (x)
is not guaranteed if just assuming f ′ (x) < 0, so, even if the original distribution is monotonic decreasing, the hl2 kernel distribution might be non-monotonic. Also, unlike the pairwise difference distribution, if the original distribution is unimodal, the pairwise mean distribution might be non-unimodal, as demonstrated by a counterexample given by Chung in 1953 and mentioned by Hodges and Lehmann in 1954 [START_REF] Hodges | Matching in paired comparisons[END_REF][START_REF] Chung | Sur les lois de probabilité unimodales[END_REF]. Theorem .9 implies that the violation of νth γ-orderliness within the hl k kernel distribution is also bounded, and the bound monotonically shrinks as k increases because the bound is in unit of the standard deviation of the hl k kernel distribution. If all hl k kernel distributions are νth γ-ordered and the distribution itself is νth γ-ordered and γ-U -ordered, then the distribution is called νth γ-U -ordered. The following theorem highlights the significance of symmetric distribution.

Theorem .17. Any symmetric distribution is νth U -ordered.

Proof. A random variable is symmetric about zero if and only if its characteristic function is real valued. Since the characteristic function of the average of k independent, identically distributed random variables is the product of the kth root of their individual characteristic functions : φ Ȳ (t) = k r=1 (φY r (t)) The succeeding theorem shows that the whl k kernel distribution is invariably a location-scale distribution if the original distribution belongs to a location-scale family with the same location and scale parameters. According to Theorem .18, the γ-weighted inequality for a right-skewed distribution can be modified as ∀0 ≤ ϵ0 1 ≤ ϵ0 2 ≤ 1 1+γ , WLM k,ϵ=1-(1-ϵ 0 1 )

1 k ,γ ≥ WLM k,ϵ=1-(1-ϵ 0 2 ) 1 k ,γ
, which holds the same rationale as the γ-weighted inequality defined in the last section. If the νth γ-orderliness is valid for the whl k kernel distribution, then all results in the last section can be directly implemented. From that, the binomial H-L mean (set the WA as BM) can be constructed (Figure 1), while its maximum breakdown point is ≈ 0.065 if ν = 3. A comparison of the biases of STM ϵ= ,n are nearly optimal with regards to concentration bounds for heavy-tailed distributions [START_REF] Devroye | Sub-gaussian mean estimators[END_REF][START_REF] Laforgue | On medians of (randomized) pairwise means in International Conference on Machine Learning[END_REF].

In 1958, Richtmyer introduced the concept of quasi-Monte Carlo simulation that utilizes low-discrepancy sequences, resulting in a significant reduction in computational expenses for large sample simulation [START_REF] Rd Richtmyer | A non-random sampling method, based on congruences, for" monte carlo" problems[END_REF]. Among various low-discrepancy sequences, Sobol sequences are often favored in quasi-Monte Carlo methods [START_REF] Im Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]. Building upon this principle, in 1991, Do and Hall extended it to bootstrap and found that the quasi-random approach resulted in lower variance compared to other bootstrap Monte Carlo procedures [START_REF] Ka Do | Quasi-random resampling for the bootstrap[END_REF]. By using a deterministic approach, the variance of mHLM k,n is much lower than that of MoM k,b= n k (SI Dataset S1), when k is small. This highlights the superiority of the median Hodges-Lehmann mean over the median of means, as it not only can provide an accurate estimate for moderate sample sizes, but also allows the use of quasi-bootstrap, where the bootstrap size can be adjusted as needed.

Methods

The robust location estimates presented in Figure 1 and SI Dataset S1 were obtained using large quasi-random samples [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] with sample size 3.686 million for the Weibull, gamma, Pareto, and lognormal distributions within specified kurtosis ranges as shown in Figure 1 to study the large sample performance. The standard errors of these estimators were computed by approximating the sampling distribution using 1000 pseudorandom samples of size n = 5184 for these distribution and the generalized Gaussian distributions with the parameter settings detailed in the SI Text.

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

1 a

 1 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 11, 2024 As the most fundamental problem in statistics, robust location estimation has many prominent solutions, such as the trimmed mean, Winsorized mean, Hodges-Lehmann estimator, Huber M -estimator, and median of means. Recent studies suggest that their maximum biases concerning the mean can be quite different, but the underlying mechanisms largely remain unclear. This study exploited a semiparametric method to classify distributions by the asymptotic orderliness of quantile combinations with varying breakdown points, showing their interrelations and connections to parametric distributions. Further deductions explain why the Winsorized mean typically has smaller biases compared to the trimmed mean; two sequences of semiparametric robust mean estimators emerge, particularly highlighting the superiority of the median Hodges-Lehmann mean. semiparametric | mean-median-mode inequality | asymptotic | unimodal | Hodges-Lehmann estimator

4 -median difference SQA ϵ= 1 4 -

 44 proposed a measure of skewness based on the SQA ϵ= 1 m (41). Groeneveld and Meeden (1984) (42) generalized these measures of skewness based on van Zwet's convex transformation[START_REF] Wr Van Zwet | Convex Transformations of Random Variables: Nebst Stellingen[END_REF] while exploring their properties. A distribution is called monotonically right-skewed if and only

6 ,

 6 this simplified inequality is valid within the range of1 6 < ϵ ≤ 1 1+γ , when 0 < γ ≤ 1. Then, it validates ∂ sup QA ∂ϵ ≤ 0 for the same range of ϵ and γ. The first and second formulae, when ϵ = 1 follows that sup QA(ϵ, γ) is contin- uous over [0,1 1+γ ]. Hence, ∂ sup QA ∂ϵ ≤ 0 holds for the entire range 0 ≤ ϵ ≤ 1 1+γ , when 0 ≤ γ ≤ 1, which leads to the assertion of this theorem.Let P kΥ denote the set of all continuous distributions whose moments, from the first to the kth, are all finite. For a right-skewed distribution, it suffices to consider the upper D R A F T bound. The monotonicity of sup P ∈P 2 Υ

  du, the proof of the first assertion is complete. The second assertion is established in Theorem 0.3. in the SI Text.

  which are the left and right sides of the desired inequality. Given that the Li PNAS | February 11, 2024 | vol. XXX | no. XX | 7

  symmetric. The conclusion follows immediately from the definition of νth U -orderliness and Theorem .2, .3, and .4.

Theorem . 18 .+

 18 whl k (x1 = λx1 + µ, . . . , x k = λx k + µ) = λwhl k (x1, . . . , x k ) + µ. Proof. whl k (x1 = λx1 + µ, ••• , x k = λx k + µ) µ = λwhl k (x1, ••• , x k ) + µ.

1 8, SWM ϵ= 1 8 , BWM ϵ= 1 8 , BM ν=2,ϵ= 1 8 , BM ν=3,ϵ= 1 8 , SQM ϵ= 1 8 , THLM k=2,ϵ= 1 8 , 8 ( 1 ,

 188888881 WiHLM k=2,ϵ= 1 Winsorized H-L mean), SQHLM k= 2 ln(2)-ln(3) 3 ln(2)-ln(7) SI Dataset S1), given their same breakdown points, with mHLM k= ln(2) 3 ln(2)-ln(7) ,ϵ= 1 8 exhibiting the smallest biases. Another comparison among the H-L estimator, the trimmed mean, and the Winsorized mean, all with the same breakdown point, yields the same result that the H-L estimator has the smallest biases (SI Dataset S1). This aligns with Devroye et al. (2016) and Laforgue, Clemencon, and Bertail (2019)'s seminal works that MoM k,b= n k and MoRM k,b

ε= 1 8 )Fig. 1 .

 81 Fig. 1. Standardized biases (with respect to µ) of fifteen robust location estimates (including two parametric estimators from REDS III for better comparison) on large quasi-random samples in four two-parameter right skewed unimodal distributions, as a function of the kurtosis. The methods are described in the SI Text.

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

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

| Li

| Li

D R A F T