Detecting Extreme Values with Order Statistics in Samples from Continuous Distributions

: In the subject of statistics for engineering, physics, computer science, chemistry, and earth sciences, one of the sampling challenges is the accuracy, or, in other words, how representative the sample is of the population from which it was drawn. A series of statistics were developed to measure the departure between the population (theoretical) and the sample (observed) distributions. Another connected issue is the presence of extreme values—possible observations that may have been wrongly collected—which do not belong to the population selected for study. By subjecting those two issues to study, we hereby propose a new statistic for assessing the quality of sampling intended to be used for any continuous distribution. Depending on the sample size, the proposed statistic is operational for known distributions (with a known probability density function) and provides the risk of being in error while assuming that a certain sample has been drawn from a population. A strategy for sample analysis, by analyzing the information about quality of the sampling provided by the order statistics in use, is proposed. A case study was conducted assessing the quality of sampling for ten cases, the latter being used to provide a pattern analysis of the statistics.


Introduction
Under the assumption that a sample of size n, was drawn from a certain population (x 1 , ..., x n ∈ X) with a known distribution (with known probability density function, PDF) but with unknown parameters (in number of m, {π 1 , ..., π m }), there are alternatives available in order to assess the quality of sampling.
One category of alternatives sees the sample as a whole-and in this case, a series of statistics was developed to measure the agreement between a theoretical (in the population) and observed (of the sample) distribution. This approach is actually a reversed engineering of the sampling distribution, providing a likelihood for observing the sample as drawn from the population. To do this for any continuous distribution, the problem is translated into the probability space by the use of a cumulative distribution function (CDF).
Formally, if PDF(x; (π j ) 1≤j≤m ) takes values on a domain D, then CDF is defined by Equation (1) and {p 1 , ..., p n } defined by Equation (2) is the series of cumulative probabilities associated with the drawings from the sample.
Recent uses of those statistics include [11] (CM), [12] (WU), [13] (KS), [14] (AD), and [15] (H1). Any of the above given test statistics are to be used, providing a risk of being in error for the assumption (or a likelihood to observe) that the sample ({x 1 , ..., x n }) was drawn from the population (X). Usually this risk of being in error is obtained from Monte Carlo simulations (see [16]) applied on the statistic in question and, in some of the fortunate cases, there is also a closed-form expression (or at least, an analytic expression) for CDF of the statistic available as well. In the less fortunate cases, only 'critical values' (values of the statistic for certain risks of being in error) for the statistic are available.
The other alternative in assessing the quality of sampling refers to an individual observation in the sample, specifically the less likely one (having associated q 1 or q n with the notations given in Equation (4)). The test statistic is g1 [15], given in Equation (11).
It should be noted that 'taken as a whole' refers to the way in which the information contained in the sample is processed in order to provide the outcome. In this scenario ('as a whole'), the entirety of the information contained in the sample is used. As it can be observed in Equations (5)-(10), each formula uses all values of sorted probabilities ({q 1 , ..., q n }) associated with the values ({x 1 , ..., x n }) contained in the sample, while, as it can be observed in Equation (11), only the extreme value (max({q 1 , ..., q n }) or min({q 1 , ..., q n })) is used; therefore, one may say that only an individual observation (the extremum portion of the sample) yields the statistical outcome.
The statistic defined by Equation (11) no longer requires cumulative probabilities to be sorted; one only needs to find the most departed probability from 0.5-see Equation (11)-or, alternatively, to find the smallest (one having associated q 1 defined by Equation (4)) and the largest (one having associated q n defined by Equation (4)), and to find which deviates from 0.5 the most (g1 Statistic = max{|q 1 − 0.5|, |q n − 0.5|}).
We hereby propose a hybrid alternative, a test statistic (let us call it TS) intended to be used in assessing the quality of sampling for the sample, which is mainly based on the less likely observation in the sample, Equation (12).
The aim of this paper is to characterize the newly proposed test statistic (TS) and to analyze its peculiarities. Unlike the test statistics assessing the quality of sampling for the sample taken as a whole (Equations (5)- (10), and like the test statistic assessing the quality of sampling based on the less likely observation of the sample, Equation (11), the proposed statistic, Equation (12), does not require that the values or their associated probabilities ({p 1 , ..., p n }) be sorted (as {q 1 , ..., q n }); since (like the g1 statistic) it uses the extreme value from the sample, one can still consider it a sort of OS [17]. When dealing with extreme values, the newly proposed statistic, Equation (12), is a much more natural construction of a statistic than the ones previously reported in the literature, Equations (5)-(10), since its value is fed mainly from the extreme value in the sample (see the max function in Equation (12)). Later, it will be given a pattern analysis, revealing that it belongs to a distinct group of statistics that are more sensitive to the presence of extreme values. A strategy of using the pool of OS (Equations (5)-(12)) including TS in the context of dealing with extreme values is given, and the probability patterns provided by the statistics are analyzed.
The rest of the paper is organized as follows. The general strategy of sampling a CDF from an OS and the method of combining probabilities from independent tests are given in Section 2, while the analytical formula for the proposed statistic is given in Section 3.1, and computation issues and proof of fact results are given in Section 3.2. Its approximation with other functions is given in Section 3.3. Combining its calculated risk of being in error with the risks from other statistics is given in Section 3.4, while discussion of the results is continued with a cluster analysis in Section 3.5, and in connection with other approaches in Section 3.6. The paper also includes an appendix of the source codes for two programs and accompanying Supplementary Material.

Addressing the Computation of CDF for OS(s)
A method of constructing the observed distribution of the g1 statistic, Equation (11), has already been reported elsewhere [15]. A method of constructing the observed distribution of the Anderson-Darling (AD) statistic, Equation (9), has already been reported elsewhere [17]; the method for constructing the observed distribution of any OS via Monte Carlo (MC) simulation, Equations (5)- (12), is described here and it is used for TS, Equation (12).
Let us take a sample size of n. The MC simulation needs to generate a large number of samples (let the number of samples be m) drawn from uniform continuous distribution ({p 1 , ..., p n } in Equation (2)). To ensure a good quality MC simulation, simply using a random number generator is not good enough. The next step (Equations (10)- (12) do not require this) is to sort the probabilities to arrive at {q 1 , ..., q n } from Equation (4) and to calculate an OS (an order statistic) associated with each sample. Finally, this series of sample statistics ({OS 1 , ..., OS w } in Figure 1) must be sorted in order to arrive at the population emulated distribution. Then, a series of evenly spaced points (from 0 to 1000 in Figure 1) corresponding to fixed probabilities (from InvCDF 0 = 0 to InvCDF 1000 = 1 in Figure 1) is to be used saving the (OS statistic, its observed CDF probability) pairs ( Figure 1).
Step 1 Step 2 Step 3 Step 4 The main idea is how to generate a good pool of random samples from a uniform U(0, 1) distribution. Imagine a (pseudo) random number generator, Rand, is available, which generates numbers from a uniform U(0, 1) distribution, from a [0, 1) interval; such an engine is available in many types of software and in most cases, it is based on Mersenne Twister [18]. What if we have to extract a sample of size n = 2? If we split in two the [0, 1) interval (then into [0, 0.5) and [0.5, 1)) then for two values (let us say v1 and v2), the contingency of the cases is illustrated in Figure 2.
An even better alternative is to do only 3 (=2 + 1) drawings (v1 + v2, 'undistinguishable'), which is ("0") to extract both from [0, 0.5), then "1") to extract one (let us say first) from [0, 0.5), and another (let us say second) from [0.5, 1), and finally, ("2") to extract both from [0.5, 1) and to keep a record for their occurrences (1, 2, 1), as well. For n numbers (Figure 3), it can be from [0, 0.5) from 0 to n of them, with their occurrences being accounted for. According to the formula given in Figure 3, for n numbers to be drawn from [0, 1), a multiple of n + 1 drawings must be made in order to maintain the uniformity of distribution (w from Figure 1 becomes n + 1). In each of those drawings, we actually only pick one of n (random) numbers (from the [0, 1) interval) as independent. In the (j + 1)-th drawing, the first j of them are to be from [0, 0.5), while the rest are to be from [0.5, 1). The algorithm implementing this strategy is given as Algorithm 1.
Algorithm 1 is ready to be used to calculate any OS (including the TS first reported here). For each sample drawn from the U(0, 1) distribution (the array v in Algorithm 1), the output of it (the array u and its associated frequencies n!/j!/(n − j)!) can be modified to produce less information and operations (Algorithm 2). Calculation of the OS (OSj output value in Algorithm 2) can be made to any precision, but for storing the result, a single data type (4 bytes) is enough (providing seven significant digits as the precision of the observed CDF of the OS). Along with a byte data type (j output value in Algorithm 2) to store each sampled OS, 5 bytes of memory is required, and the calculation of n!/(n − j)!/j! can be made at a later time, or can be tabulated in a separate array, ready to be used at a later time.
Output OSj, j EndFor Output data: (n+1) OS and their occurrences As given in Algorithm 2, each use of the algorithm sampling OS will produce two associated arrays: OSj (single data type) and j (byte data type); each of them with n + 1 values. Running the algorithm r0 times will require 5 · (n + 1) · r0 bytes for storage of the results and will produce (n + 1) · r0 OSs, ready to be sorted (see Figure 1). With a large amount of internal memory (such as 64 GB when running on a 16/24 cores 64 bit computers), a single process can dynamically address very large arrays and thus can provide a good quality, sampled OS. To do this, some implementation tricks are needed (see Table 1).

Constant/Variable/Type Value Meaning
stt ← record v:single; c:byte; end (OSj, j) pair from Algorithm 2 stored in 5 bytes mem ← 12,800,000,000 in bytes, 5*mem ← 64Gb, hardware limit buf ← 1,000,000 the size of a static buffer of data (5*buf bytes) stst ← array[0..buf-1]of stt static buffer of data dyst ← array of stst dynamic array of buffers lvl ← 1000 lvl + 1: number of points in the grid (see Figure 1) Depending on the value of the sample size (n), the number of repetitions (r2) for sampling of OS, using Algorithm 2, from r0 ← mem/(n + 1) runs, is r2 ← r0 · (n + 1), while the length (sts) of the variable (CDFst) storing the dynamic array (dyst) from Table 1 is sts ← 1 + r2/bu f . After sorting the OSs (of sttype, see Table 1; total number of r2) another trick is to extract a sample series at evenly spaced probabilities from it (from InvCDF 0 to InvCDF 1000 in Figure 1). For each pair in the sample (lvli varying from 0 to lvl = 1000 in Table 1), a value of the OS is extracted from CDFst array (which contains ordered OS values and frequencies indexed from 0 to r2−1), while the MC-simulated population size is r0 · 2 n . A program implementing this strategy is available upon request (project_OS.pas).
The associated objective (with any statistic) is to obtain its CDF and thus, by evaluating the CDF for the statistical value obtained from the sample, Equations (5)- (12), to associate a likelihood for the sampling. Please note that only in the lucky cases is it possible to do this; in the general case, only critical values (values corresponding to certain risks of being in error) or approximation formulas are available (see for instance [1][2][3]5,[7][8][9]). When a closed form or an approximation formula is assessed against the observed values from an MC simulation (such as the one given in Table 1), a measure of the departure such as the standard error (SE) indicates the degree of agreement between the two. If a series of evenly spaced points (lvl + 1 points indexed from 0 to lvl in Table 1) is used, then a standard error of the agreement for inner points of it (from 1 to lvl − 1, see Equation (13)) is safe to be computed (where p i stands for the observed probability whilep i for the estimated one).
In the case of lvl + 1, evenly spaced points in the interval [0, 1] in the context of MC simulation (as the one given in Table 1) providing the values of OS statistic in those points (see Figure 1), the observed cumulative probability should (and is) taken as p i = i/lvl, whilep i is to be (and were) taken from any closed form or approximation formula for the CDF statistic (labeledp) aŝ p i =p(InvCDF i ), where InvCDF i are the values collected by the strategy given in Figure 1 operating on the values provided by Algorithm 2. Before giving a closed form for CDF of TS (Equation (12)) and proposing approximation formulas, other theoretical considerations are needed.

Further Theoretical Considerations Required for the Study
When the PDF is known, it does not necessarily imply that its statistical parameters ((π j ) 1≤j≤m in Equations (1)-(3)) are known, and here, a complex problem of estimating the parameters of the population distribution from the sample (it then uses the same information as the one used to assess the quality of sampling) or from something else (and then it does not use the same information as the one used to assess the quality of sampling) can be (re)opened, but this matter is outside the scope of this paper.
The estimation of distribution parameters (π j ) 1≤j≤m for the data is, generally, biased by the presence of extreme values in the data, and thus, identifying the outliers along with the estimation of parameters for the distribution is a difficult task operating on two statistical hypotheses. Under this state of facts, the use of a hybrid statistic, such as the proposed one in Equation (12), seems justified. However, since the practical use of the proposed statistics almost always requires estimation of the population parameters (and in the examples given below, as well), a certain perspective on estimation methods is required.
Assuming that the parameters are obtained using the maximum likelihood estimation method (MLE, Equation (14); see [19]), one could say that the uncertainty accompanying this estimation is propagated to the process of detecting the outliers. With a series of τ statistics (τ = 6 for Equations (5)-(10) and τ = 8 for Equations (5)-(12)) assessing independently the risk of being in error (let be α 1 , ..., α τ those risks), assuming that the sample was drawn from the population, the unlikeliness of the event (α FCS in Equation (15) below) can be ascertained safely by using a modified form of Fisher's "combining probability from independent tests" method (FCS, see [10,20,21]; Equation (15) Two known symmetrical distributions were used (PDF, see Equation (1)) to express the relative deviation from the observed distribution: Gauss (G2 in Equation (16)) and generalized Gauss-Laplace (GL in Equation (17)), where (in both Equations (16) and (17) The distributions given in Equations (16) and (17) will be later used to approximate the CDF of TS as well as in the case studies of using the order statistics. For a sum (x ← p 1 +...+ p n in Equation (18)) of uniformly distributed (p 1 , ..., p n ∈ U(0, 1)) deviates (as {p 1 , ..., p n } in Equation (2)) the literature reports the Irwin-Hall distribution [22,23]. The CDF IH (x; n) is:

The Analytical Formula of CDF for TS
The CDF of TS depends (only) on the sample size (n), e.g., CDF TS (x; n). As the proposed equation, Equation (12), resembles (as an inverse of) a sum of normal deviates, we expected that the CDF TS will also be connected with the Irwin-Hall distribution, Equation (18). Indeed, the conducted study has shown that the inverse (y ← 1/x) of the variable (x) following the TS follows a distribution (1/TS) of which the CDF is given in Equation (19). Please note that the similarity between Equations (18) and (19) is not totally coincidental; 1/TS (see Equation (12)) is more or less a sum of uniform distributed deviates divided by the highest one. Also, for any positive arbitrary generated series, its ascending (x) and descending (1/x) sorts are complementary. With the proper substitution, CDF 1/TS (y; n) can be expressed as a function of CDF IH -see Equation (20).
Unfortunately, the formulas, Equation (18) to Equation (20), are not appropriate for large n and p (p = CDF 1/TS (y; n) from Equation (19)), due to the error propagated from a large number of numerical operations (see further Table 2 in Section 3.2). Therefore, for p > 0.5, a similar expression providing the value for α = 1 − p is more suitable. It is possible to use a closed analytical formula for α = 1 − CDF 1/TS (y; n) as well, Equation (21). Equation (21) resembles the Irwin-Hall distribution even more closely than Equation (20)-see Equation (22).
For consistency in the following notations, one should remember the definition of CDF, see Equation (1), and then we mark the connection between notations in terms of the analytical expressions of the functions, Equation (23): since InvCDF TS (p; n)·InvCDF 1/TS (p; n) = 1.
One should notice (Equation (1); Equation (23)) that the infimum for the domain of 1/TS (1) is the supremum for the domain of TS (1) and the supremum (n) for the domain of 1/TS is the infimum (1/n) for the domain of TS. Also, TS has the median (p = α = 0.5) at 2/(n + 1), while 1/TS has the median (which is also the mean and mode) at (n + 1)/2. The distribution of 1/TS is symmetrical.
For n = 2, the p = CDF 1/TS (y; n) is linear (y + p = 2), while for n = 3, it is a mixture of two square functions: 2p = (3 − y) 2 , for p ≤ 0.5 (and y ≥ 2), and 2p + (y − 1) 2 = 1 for p ≥ 0.5 (and x ≤ 2). With the increase of n, the number of mixed polynomials of increasing degree defining its expression increases. Therefore, it has no way to provide an analytical expression for InvCDF of 1/TS, not even for certain p values (such as 'critical' analytical functions).
The distribution of 1/TS can be further characterized by its central moments (Mean µ, Variance σ 2 , Skewness γ 1 , and Kurtosis κ in Equation (24)), which are closely connected with the Irwin-Hall distribution.

Computations for the CDF of TS and Its Analytical Formula
Before we proceed in providing the simulation results, some computational issues must be addressed. Any of the formulas provided for CDF of TS (Equations (19) and (21); or Equations (20) and (22) both connected with Equation (18)), will provide almost exact calculations as long as computations with the formulas are conducted with an engine or package that performs the operations with rational numbers to an infinite precision (such as is available in the Mathematica software [24]), when also the value of y (y ← 1/x, of floating point type) is converted to a rounded, rational number. Otherwise, with increasing n, the evaluation of CDF for TS using either Equation (19) to Equation (22) carries huge computational errors (see the alternating sign of the terms in the sums of Equations (18), (19), and (21)). In order to account for those computational errors (and to reduce their magnitude) an alternate formula for the CDF of TS is proposed (Algorithm 3), combining the formulas from Equations (19) and (21), and reducing the number of summed terms.

Algorithm 3: Avoiding computational errors for TS.
Input data: n (n ≥ 2, integer), x (1 ≤ x ≤ 1/n, real number, double precision) y ← 1/x; //p 1/TS ← Equation (19), Output data: α = α 1/TS =p TS ← CDF TS (x; n) and p = p 1/TS =α TS ← 1− p TS Table 2 contains the sums of the residuals (SS = ∑ 999 i=1 (p i −p i ) 2 in Equation (13), lvl = 1000) of the agreement between the observed CDF of TS (p i = i/1000, for i from 1 to 999) and the calculated CDF of TS (thep i values are calculated using Algorithm 3 from x i = InvCDF(i/1000; n) for i from 1 to 999) for some values of the sample size (n). To prove the previous given statements, Table 2 provides the square sums of residuals computed using three alternate formulas (from Equation (20) and from Equation (22), along with the ones from Algorithm 3). In red: computing affected digits.
As given in Table 2, the computational errors by using either Equation (20) (or Equation (19)) and Equation (22) (or Equation (21)) until n = 34 are reasonably low, while from n = 42, they become significant. As can be seen (red values in Table 2), double precision alone cannot cope with the large number of computations, especially as the terms in the sums are constantly changing their signs (see (−1) k in Equations (19) and (21)).
The computational errors using Algorithm 3 are reasonably low for the whole domain of the simulated CDF of TS (with n from 2 to 55), but the combined formula (Algorithm 3) is expected to lose its precision for large n values, and therefore, a solution to safely compute (CDF for I H, TS and 1/TS) is to operate with rational numbers.
One other alternative is to use GNU GMP (Multiple Precision Arithmetic Library [25]). The calculations are the same (Algorithm 3); the only difference is the way in which the temporary variables are declared (instead of double, the variables become mp f _t initialized later with a desired precision).
Regarding Table 2, Algorithm 4 listed data, from n = 2 to n = 55, the calculation of the residuals were made with double (64 bits), extended (FreePascal 80 bits), and mp f _t-128 bits (GNU GMP). The sum of residuals (for all n from 2 to 55) differs from double to extended with less than 10 −11 and the same for mp f_t with 128 bits, which safely provides confidence in the results provided in Table 2 for the combined formula (last column, Algorithm 4). The deviates for agreement in the calculation of CDF for TS are statistically characterized by SE (Equation (13)), min, and max in Table 3.

Approximations of CDF of TS with Known Functions
Considering, once again, Equation (24), for sufficiently large n, the distribution of 1/TS is approximately normal (Equation (26). For normal Gauss distribution, see Equation (16)). PDF 1/TS (y; n) n→∞ −−→PDF G2 ((n+1)/2; (n−1)/12)). (26) Even better (than Equation (26)), for large values of n, a generalized Gauss-Laplace distribution (see Equation (17)) can be used to approximate the 1/TS statistic. Furthermore, for those looking for critical values of the TS statistic, the approximation of the 1/TS statistic to a generalized Gauss-Laplace distribution may provide safe critical values for large n. One way to derive the parameters of the generalized Gauss-Laplace distribution approximating the 1/TS statistic is by connecting the kurtosis and skewness of the two (Equation (27)).
As can be observed in Figure 9, the confidence in approximation of 1/TS with the GL increases with the sample size (n), but the increase is less than linear. The tendency is to approximately linearly decrease with an exponential increase.  (29)).
The calculation of CDF for 1/TS is a little tricky, as anticipated previously (see Section 3.2). To avoid the computation errors in the calculation of CDF TS , a combined formula is more appropriate (Algorithms 3 and 4). With p 1/TS ← CDF 1/TS (y; n) and α 1/TS ← 1 − CDF 1/TS (y; n), depending on the value of y (y ← 1/x, where x is the sample statistic of TS, Equation (12)), only one (from α and p, where α + p = 1) is suitable for a precise calculation.

The Use of CDF for TS to Measure the Departure between an Observed Distribution and a Theoretical One
With any of Equations (5)-(12), a likelihood to observe an observed sample can be ascertained. One may ask which statistic is to be trusted. The answer is, at the same time, none and all, as the problem of fitting the data to a certain distribution involves the estimation of the distribution's parameters-such as using MLE, Equation (14). In this process of estimation, there is an intrinsic variability that cannot be ascertained by one statistic alone. This is the reason that calculating the risk of being in error from a battery of statistics is necessary, Equation (15).
Also, one may say that the g1 statistic (Equation (11)) is not associated with the sample, but to its extreme value(s), while others may say the opposite. Again, the truth is that both are right, as in certain cases, samples containing outliers are considered not appropriate for the analysis [27], and in those cases, there are exactly two modes of action: to reject the sample or to remove the outlier(s). Figure 10 gives the proposed strategy of assessing the samples using order statistics. As other authors have noted, in nonparametric problems, it is known that order statistics, i.e., the ordered set of values in a random sample from least to greatest, play a fundamental role. 'A considerable amount of new statistical inference theory can be established from order statistics assuming nothing stronger than continuity of the cumulative distribution function of the population' as [28] noted, a statement that is perfectly valid today.
In the following case studies, the values of the sample statistics were calculated with Equations (5)-(10) (AD, KS, CM, KV, WU, H1; see also Figure 10), while the risks of being in error-associated with the values of sample statistics (α Statistic for those)-were calculated with the program developed and posted online available at http://l.academicdirect.org/Statistics/tests. The g1 Statistic (Equation (11)) and α g1 were calculated as given in [15], while the TS Statistic (Equation (12)) was calculated with Algorithm 4. For FCS and α FCS , Equation (15) was used.

The Patterns in the Order Statistics
A cluster analysis on the risks of being in error, provided by the series of order statistics on the case studies considered in this study, may reveal a series of peculiarities (Figures 11 and 12). The analysis given here is based on the series of the above given case studies in order to illustrate similarities (and not to provide a 'gold standard' as in [32] or in [33]). Both clustering methods illustrated in Figures 11 and 12 reveal two distinct groups of statistics: {AD, CM, KV, WU, KS} and {H1, TS, g1}. The combined test FCS is also attracted (as expected) to the largest group. When looking at single Euclidean distances (Figure 11) of the largest group, two other associations should be noticed {AD, CM, KS} and {KV, WU}, suggesting that those groups carry similar information, but when looking at the Pearson disagreements (Figure 12), we must notice that the subgroups are changed {CM, KV, WU}, {AD}, and {KS}, with no hint of an association with their calculation formulas (Equations (5)-(9)); therefore, their independence should not be dismissed. The second group {H1, TS, g1} is more stable, maintaining the same clustering pattern of the group ({H1, TS}, {g1} in Figure 12).
Taking into account that the g1 test (Equation (11)) was specifically designed to account for outliers suggests that the H1 and TS tests are more sensitive to the outliers than other statistics, and therefore, when the outliers (or just the presence of extreme values) are the main concern in the sampling, it is strongly suggested to use those tests. The H1 statistic is a Shannon entropy formula applied in the probability space of the sample. When accounting for this aspect in the reasoning, the rassociation of the H1 with TS suggests that TS is a sort of entropic measure (max-entropy, to be more exact [34], a limit case of generalized Rényi's entropy [35]). Again, the g1 statistic is alone in this entropic group, suggesting that it carries a unique fingerprint about the sample-specifically, about its extreme value (see Equation (11))-while the others account for the context (the rest of the sampled values, Equations (10) and (12)).
Regarding the newly proposed statistic (TS), from the given case studies, the fact that it belongs to the {H1, TS, g1} group strongly suggests that it is more susceptible to the presence of outliers (such as g1, purely defined for this task, and unlike the well known statistics defined by Equations (5)- (9)).
Moreover, one may ask that, if based on the risks being in error provided by the statistics from case studies 1 to 10, some peculiarity about TS or another statistic involved in this study could be revealed. An alternative is to ask if the values of risks can be considered to be belonging to the same population or not, and for this, the K-sample Anderson-Darling test can be invoked [36]. With the series of probabilities, there are actually 2 9 − 1 − 9 = 502 tests to be conducted (for each subgroup of 2, 3, 4, 5, 6, 7, 8, and 9 statistics picked from nine possible choices) and for each of them, the answer is same: At the 5% risk of being in error, it cannot be rejected that the groups (of statistics) were selected from identical populations (of statistics), so, overall, any of those statistics perform the same.
The proposed method may find its uses in testing symmetry [37], as a homogeneity test [38] and, of course, in the process of detecting outliers [39].

Another Rank Order Statics Method and Other Approaches
The series of rank order statistics included in this study, Equations (5)- (11), covers the most known rank order statistics reported to date. However, when considering a new order statistic not included there, the use of it in the context of combining methods, Equation (15), only increases the degrees of freedom τ, while the design of using ( Figure 10) is changed accordingly.
It should be noted that the proposed approach is intended to be used for small sample sizes, when no statistic alone is capable of high precision and high trueness. With the increasing sample size, all statistics should converge to the same risk of being in error and present other alternatives, such as the superstatistical approach [40]. In the same context, each of the drawings included in the sample are supposed to be independent. In the presence of correlated data (such as correlated in time), again, other approaches, such as the one communicated in [41], are more suited.

Conclusions
A new test statistic to be used to measure the agreement between continuous theoretical distributions and samples drawn from TS was proposed. The analytical formula of the TS cumulative distribution function was obtained. The comparative study against other order statistics revealed that the newly proposed statistic carries distinct information regarding the quality of the sampling. A combined probability formula from a battery of statistics is suggested as a more accurate measure for the quality of the sampling. Therefore Equation (15) combining the probabilities (the risks of being in error) from Equation (5) to Equation (12) is recommended anytime when extreme values are suspected being outliers in samples from continuous distributions. and to provide the data for Figure 8. FreePascal (with GNU GMP, freeware) were used to assess numerically the agreement for TS statistic (Tables 2 and 3, Figure 8). StatSoft Statistica (v.7, licensed) was used to obtain Figures 11 and 12.

Conflicts of Interest:
The author declares no conflict of interest.