Gastwirth’s location estimator

The problem of outliers – data points that are substantially inconsistent with the majority of the other points in a dataset – arises frequently in the analysis of numerical data.  The practical importance of outliers lies in the fact that even a few of these points can badly distort the results of an otherwise reasonable data analysis.  This outlier-sensitivity problem is often particularly acute for classical data characterizations and analysis methods like means, standard deviations, and linear regression analysis.  As a consequence, a range of outlier-resistant methods have been developed for many different applications, and new methods continue to be developed.  For example, the R package robustbase that I have discussed in previous posts includes outlier-resistant methods for estimating location (i.e., outlier-resistant alternatives to the mean), estimating scale (outlier-resistant alternatives to the standard deviation), quantifying asymmetry (outlier-resistant alternatives to the skewness), and fitting regression models.  In Exploring Data in Engineering, the Sciences, and Medicine, I discuss a number of outlier-resistant methods for addressing some of these problems, including Gastwirth’s location estimator, an alternative to the mean that is the subject of this post.

The mean is the best-known location estimator, and it gives a useful assessment of the “typical” value of any numerical sequence that is reasonably symmetrically distributed and free of outliers.  The outlier-sensitivity of the mean is severe, however, which motivates the use of outlier-resistant alternatives like the median.  While the median is almost as well-known as the mean and extremely outlier-resistant, it can behave unexpectedly (i.e., “badly”) as a result of its non-smooth character.  This point is illustrated in Fig. 7.23 in Exploring Data, identical in character to the figure shown below (this figure is slightly different because it uses a different seed to generate the random numbers on which it is based).  Specifically, this plot shows a sequence of 201 data points, constructed as follows.  The first 100 points are normally distributed with mean 1 and standard deviation 0.1, the 101^st point is equal to zero, and points 102 through 201 are normally distributed with mean -1 and standard deviation 0.1.  Small changes in this dataset in the specific form of deleting points can result in very large changes in the computed median.  Specifically, in this example, the first 100 points lie between 0.768 and 1.185 and the last 100 points lie between -0.787 and -1.282; because the central data point lies between these two equal-sized groups, it defines the median, which is 0.  The mean is quite close to this value, at -0.004, but the situation changes dramatically if we omit either the first two or the last two points from this data sequence.  Specifically, the median value computed from points 1 through 199 is 0.768, while that computed from points 3 through 201 is -0.787.  In contrast, the mean values for these two modified sequences are 0.006 and -0.014.  Thus, although the median is much less sensitive than the mean to contamination from outliers, it is extremely sensitive to the 1% change made in this example for this particular dataset. 

The fact that the median is not “universally the best location estimator” provides a practical motivation for examining alternatives that are intermediate in behavior between the very smooth but very outlier-sensitive mean and the very outlier-insensitive but very non-smooth median.  Some of these alternatives were examined in detail in the book Robust Estimates of Location: Survey and Advances, by D.F. Andrews, P.J. Bickel, F.R. Hampel, P.J. Huber, W.H. Rogers, and J.W. Tukey, published by Princeton University Press in 1972 (according to the publisher's website, this book is out of print, but used copies are available through distributors like Amazon or Barnes and Noble).  The book summarizes the results of a year-long study of 68 different location estimators, including both the mean and the median.  The fundamental criteria for inclusion in this study were, first, that the estimators had to be computable from any given sequence of real numbers, and second, that they had to be both location and scale-invariant.  Specifically, if a given data sequence {x_k} yielded a result m, the scaled and shifted data sequence {Ax_k + b} should yield the result Am+b, for any numbers A and b.  The study was co-authored by six statistical researchers with differing opinions and points of view, but two of the authors – D.F. Andrews and F.R. Hampel – included the Gastwirth estimator (described in detail below) in their list of favorites.  For example, Hampel characterized this estimator as one of a small list of those that were “never bad at the distributions considered.”  Also, in contrast to many of the location estimators considered in the study, Gastwirth’s estimator does not require iterative computations, making it simpler to implement.

Specifically, Gastwirth’s location estimator is a weighted sum of three order statistics.  That is, to compute this estimator, we first sort the data sequence in ascending order.  Then, we take the values that are one-third of the way up this sequence (the 0.33 quantile), half way up the sequence (i.e., the median, or 0.50 quantile), and two-thirds of the way up the sequence (the 0.67 quantile).  Given these three values, we then form the weighted average, giving the central (median) value a weight of 40% and the two extreme values each a weight of 30%.  This is extremely easy to do in R, with the following code:

Gastwirth <- function(x,...){

  #

  ordstats = quantile(x, probs=c(1/3,1/2,2/3),...)

  wts = c(0.3,0.4,0.3)

  sum(wts*ordstats)

  #

}

The key part of this code is the first line, which computes the required order statistics (i.e., the quantiles 1/3, 1/2, and 2/3) using the built-in quantile function.  The first argument passed to this function is x, the vector of data values to be characterized, and the second argument (probs) defines the specific quantiles we wish to compute.  The ellipses in the Gastwirth procedure’s command line is passed to the quantile function; several parameters are possible (type “help(quantile)” in your R session for details), but one of the most useful is na.rm, a logical variable that specifies how missing data values are to be handled.  The default is “FALSE” and this causes the Gastwirth procedure to return the missing data value “NA” if any values of x are missing; the alternative “TRUE” computes the Gastwirth estimator from the non-missing values, giving a numerical result.  The three-element vector wts defines the quantile weights that define the Gastwirth estimator, which the final sum statement computes.

For the data example considered above, the Gastwirth estimator yields the location estimate -0.001 for the complete dataset, 0.308 for points 1 to 199 (vs. 0.768 for the median), and -0.317 for points 3 to 201 (vs. -0.787 for the median).  Thus, while it does not perform nearly as well as the mean for this example, it performs substantially better than the median. 

For the infinite-variance Cauchy distribution that I have discussed in several previous posts, the Gastwirth estimator performs similarly to the median, yielding a useful estimate of the center of the data distribution, in contrast to the mean, which doesn’t actually exist for this distribution (that is, the first moment does not exist for the Cauchy distribution).  Still, the distribution is symmetric about zero, so the median is well-defined, as is the Gastwirth estimator, and both should be zero for this distribution.  The above figure shows the results of applying these three estimators – the mean, the median, and Gastwirth’s estimator – to 1,000 independent random samples drawn from the Cauchy distribution.  Specifically, this figure gives a boxplot summary of these results, truncated to the range from -3 to 3 to show the range of variation of the median and Gastwirth estimator (without this restriction, the boxplot comparison would be fairly non-informative, since the mean values range from approximately -161 to 27,793, reflecting the fact that the mean is not a consistent location estimator for the Cauchy distribution).   To generate these results, the replicate function in R was used, followed by the apply function, as follows:

    RandomSampleFrame = replicate(1000, rt(n=200,df=1))
    BoxPlotVector = apply(RandomSampleFrame, MARGIN=2, Gastwirth)

The replicate function creates a data frame with the number of columns specified by the first argument (here, 1000), and each column generated by the R statement that appears as the second argument.  In this case, this second argument is the command rt, which generates a sequence of n statistically independent random numbers drawn from the Student’s t-distribution with the number of degrees of freedom specified by the df argument (here, this is 1, corresponding to the fact that the Cauchy distribution is the Student’s t-distribution with 1 degree of freedom).   Thus, RandomSampleFrame is a data frame with 200 rows and 1,000 columns, each of which may be regarded as a Cauchy-distributed random sample.  The apply function applies the function specified in the third argument (here, the Gastwirth procedure listed above) to the columns (MARGIN=2 specifies columns; MARGIN=1 would specify rows) of the data frame specified in the first argument.  The result is BoxPlotVector, a vector of 1,000 Gastwirth estimates, one for each random sample generated by the replicate function above.

At the other extreme, in the limit of infinite degrees of freedom, the Student’s t-distribution approaches a Gaussian limit.  The figure above shows the same comparison as before, except for the Gaussian distribution instead of the Cauchy distribution.  Here, the mean is the best possible location estimator and it clearly performs the best, but the point of this example is that Gastwirth’s location estimator performs better than the median.  In particular, the interquartile distance (i.e., the width of the “box” in each boxplot) for the mean is 0.094, it is 0.113 for the median, and it is 0.106 for Gastwirth’s estimator.

Another application area where very robust estimators like the median often perform poorly is that of bimodal distributions like the arc-sine distribution whose density is plotted above.  This distribution is a symmetric beta distribution, with both shape parameters equal to 0.5 (see Exploring Data, Sec. 4.5.1 for further discussion of this distribution).  Because it is symmetrically distributed on the interval from 0 to 1, the location parameter for this distribution is 0.5 and all three of the location estimators considered here yield values that are accurate on average, but with different levels of precision.  This point is shown in the figure below, which again provides boxplot comparisons for 1,000 random samples drawn from this distribution, each of length 200, for the mean, median, and Gastwirth location estimators.  As in the Gaussian case considered above, the mean performs best here, with an interquartile distance of 0.035, the median performs worst, with an interquartile distance of 0.077, and Gastwirth’s estimator is intermediate, with an interquartile distance of 0.060.

The point of this post has been to illustrate a location estimator with properties that are intermediate between those of the much better-known mean and median.  In particular, the results presented here for the Cauchy distribution show that Gastwirth’s estimator is intermediate in outlier sensitivity between the disastrously sensitive mean and the maximally insensitive median.  Similarly, the first example demonstrated that Gastwirth’s estimator is also intermediate in smoothness between the maximally smooth mean and the discontinuous median: the sensitivity of Gastwirth’s estimator to data editing in “swing-vote” examples like the one presented here is still undesirably large, but much better than that of the median.  Finally, the results presented here for the Gaussian and arc-sine distributions show that Gastwirth’s estimator is better-behaved for these distributions than the median.  Because it is extremely easy to implement in R, Gastwirth’s estimator seems worth knowing about.

Is the “Long Tail” a Useless Concept?

In response to my last post, “The Long Tail of the Pareto Distribution,” Neil Gunther had the following comment:

           

“Unfortunately, you've fallen into the trap of using the ‘long tail’ misnomer. If you think about it, it can't possibly be the length of the tail that sets distributions like Pareto and Zipf apart; even the negative exponential and Gaussian have infinitely long tails.”

He goes on to say that the relevant concept is the “width” or the “weight” of the tails that is important, and that a more appropriate characterization of these “Long Tails” would be “heavy-tailed” or “power-law” distributions.

Neil’s comment raises an important point: while the term “long tail” appears a lot in both the on-line and hard-copy literature, it is often somewhat ambiguously defined.  For example, in his book, The Long Tail, Chris Anderson offers the following description (page 10):

           

“In statistics, curves like that are called ‘long-tailed distributions’ because the tail of the curve is very long relative to the head.”

The difficulty with this description is that it is somewhat ambiguous since it says nothing about how to measure “tail length,” forcing us to adopt our own definitions.  It is clear from Neil’s comments that the definition he adopts for “tail length” is the width of the distribution’s support set.  Under this definition, the notion of a “long-tailed distribution” is of extremely limited utility: the situation is exactly as Neil describes it, with “long-tailed distributions” corresponding to any distribution with unbounded support, including both distributions like the Gaussian and gamma distribution where the mean is a reasonable characterization, and those like the Cauchy and Pareto distribution where the mean doesn’t even exist. 

The situation is analogous to that of confidence intervals, which characterize the uncertainty inherited by any characterization computed from a collection of uncertain (i.e., random) data values.  As a specific example consider the mean: the sample mean is the arithmetic average of N observed data samples, and it is generally intended as an estimate of the population mean, defined as the first moment of the data distribution.  A q% confidence interval around the sample mean is an interval that contains the population mean with probability at least q%.  These intervals can be computed in various ways for different data characterizations, but the key point here is that they are widely used in practice, with the most popular choices being the 90%, 95% and 99% confidence intervals, which necessarily become wider as this percentage q increases.  (For a more detailed discussion of confidence intervals, refer to Chapter 9 of Exploring Data in Engineering, the Sciences, and Medicine.)  We can, in principle, construct 100% confidence intervals, but this leads us directly back to Neil’s objection: the 100% confidence interval for the mean is entire support set of the distribution (e.g., for the Gaussian distribution, this 100% confidence interval is the whole real line, while for any gamma distribution, it is the set of all positive numbers).  These observations suggest the following notion of “tail length” that addresses Neil’s concern while retaining the essential idea of interest in the business literature: we can compare the “q% tail length” of different distributions for some q less than 100.

In particular, consider the case of J-shaped distributions, defined as those like the Pareto type I distribution whose distribution p(x) decays monotonically with increasing x, approaching zero as x goes to infinity.  The plot below shows two specific examples to illustrate the idea: the solid line corresponds to the (shifted) exponential distribution:

            p(x) = e^–(x-1)

for all x greater than or equal to 1 and zero otherwise, while the dotted line represents the Pareto type I distribution with location parameter k = 1 and shape parameter a = 0.5 discussed in my last post.  Initially, as x increases from 1, the exponential density is greater than the Pareto density, but for x larger than about 3.5, the opposite is true: the exponential distribution rapidly becomes much smaller, reflecting its much more rapid rate of tail decay.

For these distributions, define the q% tail length to be the distance from the minimum possible value of x (the “head” of the distribution; here, x = 1) to the point in the tail where the cumulative probability reaches q% (i.e., the value x_q where x < x_q with probability q%).  In practical terms, the q% tail length tells us how far out we have to go in the tail to account for q% of the possible cases.  In R, this value is easy to compute using the quantile function included in most families of available distribution functions.  As a specific example, for the Pareto type I distribution, the function qparetoI in the VGAM package gives us the desired quantiles for the distribution with specified values of the parameters k (designated “scale” in the qparetoI call) and a (designated “shape” in the qparetoI call).  Thus, for the case k = 1 and a = 0.5 (i.e., the dashed curve in the above plot), the “90% tail length” is given by:

> qparetoI(p=0.9,scale=1,shape=0.5)

[1] 100

>

For comparison, the corresponding shifted exponential distribution has the 90% tail length given by:

> 1 + qexp(p = 0.9)

[1] 3.302585

>

(Note that here, I added 1 to the exponential quantile to account for the shift in its domain from “all positive numbers” – the domain for the standard exponential distribution – to the shifted domain “all numbers greater than 1”.)  Since these 90% tail lengths differ by a factor of 30, they provide a sound basis for declaring the Pareto type I distribution to be “longer tailed” than the exponential distribution.

These results also provide a useful basis for assessing the influence of the decay parameter a for the Pareto distribution.  As I noted last time, two of the examples I considered did not have finite means (a = 0.5 and 1.0), and none of the four had finite variances (i.e., also a = 1.5 and 2.0), rendering moment characterizations like the mean and standard deviation fundamentally useless.  Comparing the 90% tail lengths for these distributions, however, leads to the following results:

            a = 0.5: 90% tail length = 100.000

            a = 1.0: 90% tail length =   10.000

            a = 1.5: 90% tail length =     4.642

            a = 2.0: 90% tail length =     3.162

It is clear from these results that the shape parameter a has a dramatic effect on the 90% tail length (in fact, on the q% tail length for any q less than 100).  Further, note that the 90% tail length for the Pareto type I distribution with a = 2.0 is actually a little bit shorter than that for the exponential distribution.  If we move further out into the tail, however, this situation changes.  As a specific example, suppose we compare the 98% tail lengths. For the exponential distribution, this yields the value 4.912, while for the four Pareto shape parameters we have the following results:

            a = 0.5: 98% tail length = 2,500.000

            a = 1.0: 98% tail length =      50.000

            a = 1.5: 98% tail length =      13.572

            a = 2.0: 98% tail length =        7.071

This value (i.e., the 98% tail length) seems a particularly appropriate choice to include here since in his book, The Long Tail, Chris Anderson notes that his original presentations on the topic were entitled “The 98% Rule,” reflecting the fact that he was explicitly considering how far out you had to go into the tail of a distribution of goods (e.g., the books for sale by Amazon) to account for 98% of the sales.

Since this discussion originally began with the question, “when are averages useless?” it is appropriate to note that, in contrast to the much better-known average, the “q% tail length” considered here is well-defined for any proper distribution.   As the examples discussed here demonstrate, this characterization also provides a useful basis for quantifying the “Long Tail” behavior that is of increasing interest in business applications like Internet marketing.  Thus, if we adopt this measure for any q value less than 100%, the answer to the title question of this post is, “No: The Long Tail is a useful concept.”

The downside of this minor change is that – as the results shown here illustrate – the results obtained using the q% tail length depend on the value of q we choose.  In my next post, I will explore the computational issues associated with that choice.

The Long Tail of the Pareto Distribution

In my last two posts, I have discussed cases where the mean is of little or no use as a data characterization.  One of the specific examples I discussed last time was the case of the Pareto type I distribution, for which the density is given by:

                        p(x) = ak^a/x^a+1

defined for all x > k, where k and a are numeric parameters that define the distribution.  In the example I discussed last time, I considered the case where a = 1.5, which exhibits a finite mean (specifically, the mean is 3 for this case), but an infinite variance.  As the results I presented last time demonstrated, the extreme data variability of this distribution renders the computed mean too variable to be useful.  Another reason this distribution is particularly interesting is that it exhibits essentially the same tail behavior as the discrete Zipf distribution; there, the probability that a discrete random variable x takes its i^th value is:

                        p_i = A/i^c,

where A is a normalization constant and c is a parameter that determines how slowly the tail decays.  This distribution was originally proposed to characterize the frequency of words in long documents (the Zipf-Estoup law), it was investigated further by Zipf in the mid-twentieth century in a wide range of applications (e.g., the distributions of city sizes), and it has become the subject of considerable recent attention as a model for “long-tailed” business phenomena (for a non-technical introduction to some of these business phenomena, see the book by Chris Anderson, The Long Tail).  I will discuss the Zipf distribution further in a later post, but one of the reasons for discussing the Pareto type I distribution first is that since it is a continuous distribution, the math is easier, meaning that more characterization results are available for the Pareto distribution. 

The mean of the Pareto type I distribution is:

                        Mean = ak/(a-1),

provided a > 1, and the variance of the distribution is finite only if a > 2.  Plots of the probability density defined above for this distribution are shown above, for k = 1 in all cases, and with a taking the values 0.5, 1.0, 1.5, and 2.0.  (This is essentially the same plot as Figure 4.17 in Exploring Data in Engineering, the Sciences, and Medicine, where I give a brief description of the Pareto type I distribution.)  Note that all of the cases considered here are characterized by infinite variance, while the first two (a = 0.5 and 1.0) are also characterized by infinite means.  As the results presented below emphasize, the mean represents a very poor characterization in practice for data drawn from any of these distributions, but there are alternatives, including the familiar median that I have discussed previously, along with two others that are more specific to the Pareto type I distribution: the geometric mean and the harmonic mean. 

The plot below emphasizes the point made above about the extremely limited utility of the mean as a characterization of Pareto type I data, even in cases where it is theoretically well-defined.  Specifically, this plot compares the four characterizations I discuss here – the mean (more precisely known as the “arithmetic mean” to distinguish it from the other means considered here), the median, the geometric mean, and the harmonic mean – for 1000 statistically independent Pareto type I data sequences, each of length N = 400, with parameters k = 1 and a = 2.0.  For this example, the mean is well-defined (specifically, it is equal to 2), but compared with the other data characterizations, its variability is much greater, reflecting the more serious impact of this distribution’s infinite variance on the mean than on these other data characterizations.

To give a more complete view of the extreme variability of the arithmetic mean, boxplots of 1000 statistically independent samples drawn from all four of the Pareto type I distribution examples plotted above are shown in the boxplots below.  As before, each sample is of size N = 400 and the parameter k has the value 1, but here the computed arithmetic means are shown for the parameter values a = 0.5, 1.0, 1.5, and 2.0; note the log scale used here because the range of computed means is so large.  For the first two of these examples, the population mean does not exist, so it is not surprising that the computed values span such an enormous range, but even when the mean is well-defined, the influence of the infinite variance of these cases is clearly evident.  It may be argued that infinite variance is an extreme phenomenon, but it is worth emphasizing here that for the specific “long tail” distributions popular in many applications, the decay rate is sufficiently slow for the variance – and sometimes even the mean – to be infinite, as in these examples.

As I have noted several times in previous posts, the median is much better behaved than the mean, so much so that it is well-defined for any proper distribution.  One of the advantages of the Pareto type I distribution is that the form of the density function is simple enough that the median of the distribution can be computed explicitly from the distribution parameters.  This result is given in the fabulous book by Johnson, Kotz and Balakrishnan that I have mentioned previously, which devotes an entire chapter (Chapter 20) to the Pareto family of distributions.  Specifically, the median of the Pareto type I distribution with parameters k and a is given by:

                        Median = 2^1/ak

Thus, for the four examples considered here, the median values are 4.0 (for a = 0.5), 2.0 (for a = 1.0), 1.587 (for a = 1.5), and 1.414 (for a = 2.0).  Boxplot summaries for the same 1000 random samples considered above are shown in the plot below, which also includes horizontal dotted lines at these theoretical median values for the four distributions.  The fact that these lines correspond closely with the median lines in the boxplots gives an indication that the computed median is, on average, in good agreement with the correct values it is attempting to estimate.  As in the case of the arithmetic means, the variability of these estimates decreases monotonically as a increases, corresponding to the fact that the distribution becomes generally better-behaved as the a parameter increases.

The geometric mean is an alternative characterization to the more familiar arithmetic mean, one that is well-defined for any sequence of positive numbers.  Specifically, the geometric mean of N positive numbers is defined as the N^th root of their product.  Equivalently, the geometric mean may be computed by exponentiating the arithmetic average of the log-transformed values.  In the case of the Pareto type I distribution, the utility of the geometric mean is closely related to the fact that the log transformation converts a Pareto-distributed random variable into an exponentially distributed one, a point that I will discuss further in a later post on data transformations.  (These transformations are the topic of Chapter 12 of Exploring Data, where I briefly discuss both the logarithmic transformation on which the geometric mean is based and the reciprocal transformation on which the harmonic mean is based, described next.)   The key point here is that the following simple expression is available for the geometric mean of the Pareto type I distribution (Johnson, Kotz, and Balakrishnan, page 577):

                        Geometric Mean = k exp(1/a)

For the four specific examples considered here, these geometric mean values are approximately 7.389 (for a = 0.5), 2.718 (for a = 1.0), 1.948 (for a = 1.5), and 1.649 (for a = 2.0).  The boxplots shown below summarize the range of variation seen in the computed geometric means for the same 1000 statistically independent samples considered above.  Again, the horizontal dotted lines indicate the correct values for each distribution, and it may be seen that the computed values are in good agreement, on average.  As before, the variability of these computed values decreases with increasing a values as the distribution becomes better-behaved.

The fourth characterization considered here is the harmonic mean, again appropriate to positive values, and defined as the reciprocal of the average of the reciprocal data values.  In the case of the geometric mean just discussed, the log transformation on which it is based is often useful in improving the distributional character of data values that span a wide range.  In the case of the Pareto type I distribution – and a number of others – the reciprocal transformation on which the harmonic mean is based also improves the behavior of the data distribution, but this is often not the case.  In particular, reciprocal transformations often make the character of a data distribution much worse: applied to the extremely well-behaved standard uniform distribution, it yields the Pareto type I distribution with a = 1, for which none of the integer moments exist; similarly, applied to the Gaussian distribution, the reciprocal transformation yields a result that is both infinite variance and bimodal.  (A little thought suggests that the reciprocal transformation is inappropriate for the Gaussian distribution because it is not strictly positive, but normality is a favorite working assumption, sometimes applied to the denominators of ratios, leading to a number of theoretical difficulties.  I will have more to say about that in a future post.)  For the case of the Pareto type I distribution, the reciprocal transformation converts it into the extremely well-behaved beta distribution, and the harmonic mean has the following simple expression:

            Harmonic mean = k(1 + a^-1)

For the four examples considered here, this expression yields harmonic mean values of 3 (for a = 0.5), 2 (for a = 1.0), 1.667 (for a = 1.5), and 1.5 (for a = 2.0).  Boxplot summaries of the computed harmonic means for the 1000 simulations of each case considered previously are shown below, again with dotted horizontal lines at the theoretical values for each case.  As with both the median and the geometric mean, it is clear from these plots that the computed values are correct on average, and their variability decreases with increasing values of the a parameter.

The key point of this post has been to show that, while averages are not suitable characterizations for “long tailed” phenomena that are becoming an increasing subject of interest in many different fields, useful alternatives do exist.  For the case of the Pareto type I distribution considered here, these alternatives include the popular median, along with the somewhat less well-known geometric and harmonic means.  In an upcoming post, I will examine the utility of these characterizations for the Zipf distribution.