A Quantile Estimate Based on Local Curve Fitting

TL;DR

This paper introduces an extended regression-based method for quantile estimation, especially effective in tail regions, applicable to multiple samples, with theoretical support and real data illustrations.

Contribution

It generalizes the RAQE method to multiple samples and provides a theoretical framework for its application and weighting scheme.

Findings

Effective in tail quantile estimation

Applicable to multiple homogeneous samples

Supported by theoretical analysis

Abstract

Quantile estimation is a problem presented in fields such as quality control, hydrology, and economics. There are different techniques to estimate such quantiles. Nevertheless, these techniques use an overall fit of the sample when the quantiles of interest are usually located in the tails of the distribution. Regression Approach for Quantile Estimation (RAQE) is a method based on regression techniques and the properties of the empirical distribution to address this problem. The method was first presented for the problem of capability analysis. In this paper, a generalization of the method is presented, extended to the multiple sample scenario, and data from real examples is used to illustrate the proposed approaches. In addition, theoretical framework is presented to support the extension for multiple homogeneous samples and the use of the uncertainty of the estimated probabilities as a weighting factor in the analysis.

Tables

Table 1. Table 1 : Particle Count Found per Wafer

27	16	16	34	14	13	10	43	23	13	8	22
14	8	15	20	79	9	13	6	18	13	34	16
15	11	4	9	12	9	38	7	15	7	4	31
9	7	35	7	8	15	13	5	4	4	13	7
39	61	27	11	10	18	14	3	15	14	8	12
9	13	35	11	23	11	9	11	42	12	4	4
15	9	8	10	11	25	10	19	8	19	11	13
12	37	44	12	9	11	74	12	27	43	4	6
5	15	3	24	22	10	23	16	5	40	27	16
5	12	23	5	30	19	8	9

Table 2. Table 2 : Annual Maximum Precipitation for Stations 25081 and 25078

	1963	1964	1965	1966	1967	1968	1969	1970	1971	1972
Station 25081	158	72.5	123	58.1	90	96.2	76.8	128.6	101.8	188
Station 25078	79.5	160	252.8	129.5	94.8	240	80.8	193	86	268
	1973	1974	1975	1976	1977	1978	1979	1980	1981	1982
Station 25081	62	74.5	62.5	136	94.3	82.5	56	84.5	55	75
Station 25078	111.5	109	85	87	95	53	78	60	98	87
	1983	1984	1985	1986	1987	1988	1989	1990	1991	1992
Station 25081	58.2	82.5	72.1	224.3	89	98	116.4	127.2	83.3	82.9
Station 25078	77.7	81	90	200	152	168	64	113	66	157
	1993	1994	1995	1996	1997	1998	1999	2000	2001	2002
Station 25081	81.3	103	53.9	114.1	63.4	60.4	83.6	72	85.9	73.8
Station 25078	115	103.5	57.5	91	63	202	122	120	170	77
	2003	2004	2005	2006
Station 25081	70.9	89.4	94.2	116.9
Station 25078	81	85	119	88

Equations

$$S_n(x)=\frac{1}{n}\sum _{i=1}^n{I}_{(-\infty ,x]}(X_i),-\infty <x<\infty ,$$

$$(a_{2i},b_{2i})=((X_{(i)}+X_{(i+1)})/2,i/n),i=1,2,\dots ,n-1,$$

$$(a_{2i-1},b_{2i-1})=(X_{(i)},(i-1/2)/n),i=1,2,\dots ,n.$$

$$\sum _{i=1}^{2m-1}{w}_i(b_{(i)}-g(x\mid \theta ){)}^2$$

$$\sum _{i=2n-1l+1}^{2n-1}{w}_i(b_{(i)}-h(x\mid \phi ){)}^2$$

$$S(x)=\frac{1}{n}\sum _{i=1}^n{I}_{(-\infty ,x]}(X_i)$$

$$E[S(x)]=F(x)$$

$$Var[S(x)]=\frac{1}{n}(F(x)(1-F(x)).$$

$$sup|S(x)-F(x)|\stackrel{a.s.}{\to }0,$$

$$Q_w=\sum _{i=1}^n{w}_i(b_i-{\hat{b}}_i{)}^2.$$

$${\hat{b}}_{1i}=\frac{1}{n_1}\sum _{i=1}^{n_1}I(X_i\le a)$$

$${\hat{b}}_{2i}=\frac{1}{n_2}\sum _{i=1}^{n_2}I(Y_i\le a)$$

$$\begin{array}{rl}{\hat{b}}_i=& \frac{n_1\widehat{b_{1i}}+n_2\widehat{b_{2i}}}{n_1+n_2}\end{array}$$

Full text

Keywords: Regression techniques, weighted least squares, homogeneous samples.

AMS Subject Classifications: 62G05, 62G32.

1 Introduction

Quantile estimation is a problem presented in different fields such as quality control, risk analysis, finance, and hydrology, among other fields. Sometimes it is necessary to know the value of a certain quantile of a given distribution in order to use certain statistical tools. For these purposes, various approaches for quantile estimation have been proposed in the literature based on curve-fitting techniques, resampling techniques, or transformations.

Curve fitting techniques require the first four moments of a sample in order to estimate the quantiles of interest by using some known curve or probability distribution. Some of the most used curves are the Pearson family of curves, the Gumbel distribution, and the lognormal distribution, just to mention a few. Daniel and Wood (1999) in their book present the strengths and limitations of the method of least squares when using curve fitting techniques. The main problem with this approach is that usually the quantiles of interest are located at the extreme of the distribution, and this technique requires information from the center of the distribution. The center values get a good fit while the adjustment at the tails gets compromised.

Transformation techniques use different models to transform the original distribution to a normal distribution. Some of the most used transformations are those proposed by Box and Cox (1964) and Johnson (1949). Other authors such as Bartlett (1947), Velleman and Hoaglin (1981), and Tan et al. (2004) present different transformations and conclusions of when to use each of them. Nevertheless, a major disadvantage with these techniques is that the quantiles estimates depend on the knowledge of the right transformation. Several transformations must be tried first before finding a useful one.

Resampling techniques use different subsets of data taken randomly with replacement from a sample of interest. By using different subsets, the parameter of interest, in this case a quantile of interest, can be estimated. This resampling usually uses bootstrap and jackknife methods. Additional information about these methods can be found in Davison and Hinkley (1997), Efron (1982), and Efron and Gong (1983), among others. However, resampling techniques are known to require additional computational complexity, making the approach awkward and unwieldy for industry practitioners.

All these previous methods have shown acceptable performance for different situations, however, all of the known methods require an overall fit on the available data, giving the same weights to all observations when usually the quantiles of interest are located in the tails of the distribution.The only approaches found in the literature that deals with different weights when fitting a distribution are the weighted maximum likelihood approaches. Nevertheless, they are parametric in nature and need to assume a specific probability distribution. For instance, Field and Smith (1994) uses this approach to reduce the ”pull” given by outliers when performing an estimation. Salazar-Alvarez et al. (2016) presented a regression approach to improve the estimates when dealing with extreme values. Regression Approach for Quantile Estimation, RAQE, was introduced to solve the problem of process capability estimation based on the percentile method. In their paper, the authors show a series of Monte Carlo simulations to demonstrate the performance of the method, and a case study was presented in the automotive related industry. The problem was limited to the use of one sample and the estimation of process capability indexes.

This paper aims to generalize the methodology proposed by Salazar-Alvarez et al. (2016). Since RAQE is based on regression techniques it can be easily confused with research areas like quantile regression and extremal quantile regression; nevertheless, RAQE does not use variable conditioning. It only assumes a set of independent and identically distributed observations. Moreover, the problem of the uncertainty of a probability is incorporated into the regression techniques to improve the curve estimation. In addition, an extension of the methodology is presented for its use in the case of multiple homogeneous samples.

In this paper the generalization of RAQE is presented in Section 2, for one or multiple homogeneous samples. Examples to illustrate the applicability of the method are presented as well. Section 3 provides a theoretical discussion that supports the proposed approach. Finally, in Section 4, general conclusions about the proposed method and suggested future work are provided.

2 Methods

This section presents a methodology to estimate quantiles by using local curve fitting for single and multiple homogeneous samples. Some examples to illustrate the use of the proposed methodology are presented as well.

2.1 RAQE

Let $X=\{X_{(i)}:i=1,\ldots,n\}$ be a random vector of ordered statistics from a sample of independent and identically distributed observations, where $x=\{x_{i}:i,\ldots,n\}$ corresponds to the observed values. This implies that $x_{1}<\ldots<x_{n}$ . Then, to estimate the desired extreme quantiles use the following procedure:

Augment the empirical distribution function (e.d.f)

[TABLE]

with additional points to assist in fitting a continuous function to approximate the discrete function $S_{n}(x)$. The additional points are defined as

[TABLE]

and

[TABLE]

If the quantile of interest is located on the lower tail; then, choose a family of curves $g(x\mid\theta)$ to fit to the ordered augmented set, $\{a_{i},b_{i}\}_{i=1}^{2m-1}$ for a selected value of $m<n/2$. Use a search method to find the set of parameters $\theta$ that minimizes the weighted least square criterion

[TABLE]

where $w_{i}=\frac{1}{\hat{\sigma}_{i}^{2}}$, correspond to the weights given to each of the observations. For the case of the empirical distribution the weights are estimated as $w_{i}=n/(b_{i}(1-b_{i}))$.

If the quantile of interest is located on the upper tail, then choose a family of curves $h(x\mid\varphi)$ to fit the ordered augmented set $\{a_{i},b_{i}\}_{2n-2l+1}^{2n-1}$ for a selected value of $l<n/2$. Use a search method to find the set of parameters $\varphi$ that minimizes

[TABLE] 2. 2.

Estimate the quantiles of interest using the inverse image of the augmented empirical distribution, $\hat{X}_{p}=g^{-1}(p)$ for the lower tail and/or $\hat{X}_{p}=h^{-1}(p)$ for the upper tail, where p represents the $p^{th}$ quantile of the data.

2.1.1 The Case of Multiple Homogeneous Samples

RAQE could also be used in the case of multiple homogeneous samples. Homogeneous samples are defined as multiple samples that have the same distribution but differ on their location and/or scale parameters. These types of samples are present in industry when dealing with lines that produce similar products, for instance, a bottling company that uses the same line to fill $375$ ml and $500$ ml bottles. For these scenarios, the problem of multiple samples is reduced to the problem of one sample by combining the k-homogeneous samples into one big sample by standardizing with the corresponding sample mean and sample standard deviation, centering and re-scaling whenever there is a significance difference found in location or scale. After the quantiles of interest are obtained from the combined sample, the estimation needs to be return to the original form $x_{ir}=s_{r}z_{j}+\bar{x}_{r}$, where $\bar{x}_{r}$ is the sample mean, $s_{r}$ is the sample standard deviation, $r$ stands for the individual samples $\left(r=1,2,\ldots,k\right)$, and $i$ represents the observation within samples $\left(i=1,2,\ldots,n_{k}\right)$ and $j$ stands for the standardize observation $\left(j=1,2,\ldots,\sum_{r=1}^{k}n_{r}\right)$. For clarity, the method is presented in Figure  1.

2.2 Real Data Examples

In this section, two examples are presented using real data to demonstrate the applicability of the methodology in real scenarios.

2.2.1 Estimation of Control Limits

To illustrate the method an example given by Chou et al. (1998) is shown in this section. Measurements of particle count in wafers were recorded from a piece of equipment used in the senmiconductor industry. A random sample of size $116$ is taken from the collected data, see Table 1. The authors seek to estimate the control limits of a control chart; nevertheless, the distribution of the sample is found to be skewed. Therefore, the proposed solution is to try different transformations and use the one with the higher p-value. This will be compared with the RAQE method.

For this data set, the Johnson transformation, the logarithmic transformation, and the square root transformation were evaluated. The selected transformation method was the Johnson transformation. The sample was then transformed and the upper and lower control limits were estimated with the moving range. The estimated values for the upper and lower control limits using transformations were $82.1991$ and $2.8146$. For the nonparametric approach, the quantiles of interest to estimate the lower and upper control limit are $0.00135$ and $0.99865$, which is equivalent of having a range of $6\sigma$ under normality. The size of $m$ and $l$ was set to work with the $25\%$ of the observed lower and the upper tails, respectively, and avoid the information located at the center of the distribution where a fit is not required. The distribution used in the sample for the lower tail was a quadratic function; while on the upper tail the Gumbel cumulative distribution function was used. The estimated values obtained with RAQE were $92.3982$ and $2.8022$. It can be seen that the estimated values using both methods are relatively close. This can also be seen in Figure 2. As can be seen, the fit at the tails is not as good as the one obtained with RAQE. The mean square error (MSE) for the fit on the tails using Johnson transformation was of $0.107$ and $0.0119$ for the lower and upper tail, respectively; meanwhile, using RAQE the MSE was of $0.012$ and $0.006$.

2.2.2 Estimation of Return Periods

Flood frequency estimation is an approach used in hydrology to analyze extreme weather events. Its main concern is to estimate the probability of the occurrence of a rare event in order to estimate the risk to the population and avoid death and monetary damages related to harsh weather. One of the main steps needed is to determine the frequency distribution that better fits the data under study in order to estimate the quantile related to the return period. If different data sets are found to be homogeneous in terms of shape, but they might have different location or scale parameters, samples might be combined adequately to increase the accuracy of the curve fitting estimation. This process can be also achieved using RAQE for multiple homogeneous samples.

The data set used in this example was obtained from climatological stations that were defined as homogenous by Arellano-Lara and Escalante-Sandoval (2014). The information was obtained by the authors using the software Rapid Extractor of Climatological Information version $3$, but it could also be obtained from google earth by installing a software application from the Comisión Nacional del Agua (CONAGUA) web site. The example uses records from the annual maximum precipitation from a region located in Sinaloa, at the northwest of Mexico. The two climatological stations are $25081$, located in the city of Culiacan, and $25078$, located in the city of Rosario. The period of time selected was from $1963$ to $2006$. The data from the two stations are presented in Table 2.

The data points are arranged from oldest to newest, from left to right; therefore, the first value corresponds to the year of 1963 while the last value represents the year 2006. In addition, tests were conducted using GNU R in order to assure that the data was homogeneous. First, the data was mapped using a scatter plot, Figure 3 and a correlation test was performed, the p-value obtained was of $0.0031$, indicating that the data is significantly correlated. However, as seen in the theoretical framework in Section 3, practitioners should not be worried about putting together correlated data to perform the RAQE. Then, a paired-t test and a Levene test were conducted and the results indicated that they have different means and variances. In addition, confidence intervals using a boostrap with $1,000$ replicas were estimated for the skewness and the kurtosis. The results show that there is sufficient statistical evidence to state that their skewness and kurtosis are not significantly different.

After all the assumptions were tested, RAQE was estimated with m and l set on $25\%$ of the data. The selected model to estimate the return period was a Gumbel function, $exp\left[-exp\left(-\left(x_{i}-\theta_{1}\right)/\theta_{2}\right)\right]$, this can be seen graphically in Figure 4. The return periods estimated were of $1,000$, $100$, and $20$ years, which correspond to the quantiles $0.999$, $0.99$, and $0.95$. The estimated quantiles for the station $25081$ were $295.031,218.54,164.51$; and for the station $25078$ the quantiles were $429.51,311.14,227.51$.

3 Theoretical Framework

The quantile estimation based on local curve fitting uses the theory of the empirical distribution. A modification of the empirical distribution function, the augmented empirical distribution, was proposed to be used with existing regression techniques to find the quantiles of interest by using its inverse image. These results were already presented by Salazar-Alvarez et al. (2016) in their previous work. Nevertheless, for completeness, each of these subjects is described next in this section together with the new theoretical developments. In addition, a discussion about correlation between samples is also presented.

3.1 The Empirical Distribution

Definition 1. Let $X_{1},\ldots,X_{n}$ be a series of independent and identically distributed copies of a random variable $X$ with continuous distribution function $F(x)=P(X\leq x)$. The statistic**

[TABLE]

is defined as the empirical distribution function (e.d.f), for all $x$, $-\infty<x<\infty$.

$I_{(-\infty,x]}(\cdot)$ is an indicator function that takes values of $1$ when evaluated in $(-\infty,x]$, and zero otherwise. An important properties of the empirical distribution is that $I_{(-\infty,x]}(X_{i})$ is a Bernoulli random variable, hence $S(x)$ is the MLE of the cumulative probability $F(x)$. In addition, at any given value of $x$

[TABLE]

[TABLE]

Also, the Gilvenko-Cantelli theorem,

[TABLE]

states that $S(x)$ converges uniform to $F(x)$ as $n\rightarrow\infty$.

As demonstrated by Salazar-Alvarez et al. (2016) the augmented empirical distribution estimates points from the true continuous cumulative probability distribution function. Hence, this recollection of points can be used to fit a curve, with a least squares procedure, to find a quantile of interest. The obtained estimates are unbiased. However, Salazar-Alvarez et al. (2016) did not considered that different $b_{i}$, $i=1,\dots,n$, are heteroscedastic where the variance can be estimated as $\hat{\sigma}_{1}^{2}=b_{i}(1-b_{i})/n$. This can be addressed by adding weights of $w_{i}=1/\hat{\sigma}_{i}^{2}$ to the least squares criterion,

[TABLE]

In addition, for the case of correlated data when using multiple homogeneous samples, suppose we have two sets of independent and identically distributed random variables $X_{1},\ldots,X_{n_{1}}$ and $Y_{1},\ldots,Y_{n_{2}}$. Suppose we analyze each sample separately and the probability at a threshold $a$ is of interest. Then, for sample $1$

[TABLE]

Similarly, for sample $2$

[TABLE]

In both cases, an unbiased estimator of $F(x)$ is obtained. If the two samples are now combined in order to obtained more information about the process,then,

[TABLE]

Since the probability is a linear combination of the individual estimators, then, it still is an unbiased estimator with variance equal to $\frac{1}{n_{1}+n_{2}}\theta(1-\theta)+\frac{2n_{1}n_{2}}{\left(n_{1}+n_{2}\right)^{2}}Cov(\hat{b}_{1i},\hat{b}_{2i})$, where $\theta=P(X\leq a)$. Therefore, practitioners should feel comfortable even though the samples are not independent.

4 Conclusions

In the literature, several techniques exist to estimate quantiles by using an overall fit of the data available. However, to estimate extreme quantiles the best method would be the one that provides the best estimation at the tails of the distributions. RAQE is a new method that was presented to solve the problem of process capability analysis when the data does not represent a normal distribution. This method is now generalized so it can be used as a technique for quantile estimation.

RAQE was presented for its use with one sample or multiple homogeneous samples. As well, real case examples were presented for the problem of estimation of control limits in quality control and the estimation of return periods in hydrology. In both cases, it is noticeable that the results are similar to the results obtained with other techniques. This technique provides practitioners with a new alternative to estimate extreme quantiles. In some cases, the estimation even requires fewer steps to obtain the quantile of interest.

As part of the future work, further analysis of this method can be done using Monte Carlo simulations; nevertheless, special considerations need to be made when capturing the behavior of the tails of different samples by using only a model. In practice, practitioners have a particular sample of interest and a particular curve can be proposed by examining the data set, and they will choose the approach that best fits their distribution using popular selection approaches such as the mean absolute deviation, BIC or AIC. This situation is similar to the one presented in this work.

4.1 Acknowledgements

4.2 Funding