Fitting the tails of the distribution

Question

I was wondering if anyone knows an R function that would make a "best fit" distribution (say normal) match the data only in the tails?

Answer 1

Although usually it's not a great idea to fit a distribution only in the tails, you can construct an estimator that will do so to varying degrees. I've outlined an approach, assuming Normal data, using the Anderson-Darling statistic, and included some notes on how to extend it in various ways.

The general idea is to construct a function of the data and parameters, in this case the Anderson-Darling (AD) statistic assuming known mean and std. deviation, and then optimize it on the data set you have. The AD statistic places a lot of weight in the tails of the distribution; saving you the effort of clicking on the link, it can be written as:

$A = n \int_{-\infty}^{\infty}[F_n(x)-F(x)]^2w(x)dF(x)$

where $F_n(x)$ is the empirical cumulative distribution function, $F(x)$ is the hypothesized distribution function (where the parameters go) and $w(x)$ is a weight function:

$w(x) = 1 / [F(x)(1-F(x))]$

As you can see, differences between the empirical and hypothesized CDFs have much, much more influence on the statistic in the tails of $F$ than in the middle, relative to the influence implicit in an equal weighting:

cdf.z <- seq(0.005,0.995,by=0.01)
ad.weights <- 1/(cdf.z*(1-cdf.z))
plot(ad.weights~qnorm(cdf.z), pch=16, type="b", lwd=2,
     xlab="X ~ std. Normal", ylab="Anderson-Darling weights")

You can, of course, substitute your own weight function, e.g., a step function which weights the middle 50% of the data "0" and the rest "1", to achieve the weighting behavior you like, but some work woutld be required to construct the code for the new weight function.

The code below minimizes the AD statistic by selecting the mean and log(std. deviation); the latter avoids minor difficulties caused by the bound $\sigma > 0$.

require(MASS)
require(nortest) # only for validation of ad.statistic code

# Anderson-Darling statistic; parms[1] = mean, parms[2] = log(std. dev.)
ad.statistic <- function(parms, x) {
  y <- sort((x-parms[1])/exp(parms[2]))
  cdf.y <- pnorm(y)  
  n <- length(x)
  i <- 1:n
  stat <- (2*i-1)*log(cdf.y) + (2*(n-i)+1)*log(1-cdf.y)
  -n - mean(stat)
}

# Comparison with ad.test (validation of function)
x <- rnorm(100)
ad.statistic(c(mean(x),log(sd(x))),x)
[1] 0.6334132
ad.test(x)
 ... removing some unnecessary lines of output ...
A = 0.6334, p-value = 0.09606

# Parameter estimation and output
foo <- optim(c(mean(x),log(sd(x))), ad.statistic, method="BFGS", x=x)

cat("AD:  ",foo$par[1],"  ",exp(foo$par[2]),"\nMLE: ",mean(x),"  ",sd(x),"\n")
AD:   -0.1817909    1.145518 
MLE:  -0.2015397    1.091143 

The estimated density functions are, not surprisingly, quite similar:

Note that as your sample size goes to $\infty$, if the data really does come from the assumed distribution, the two estimators will converge (both to the true values and each other of course.) Thus with larger sample sizes the value of fitting with the tails weighted heavily would seem to be low. On the other hand, as Procrastinator has pointed out, with a small sample you're probably not going to get good results, since there's not much data in the tails.

If you wish to use a distribution other than the Normal, the only line of code that refers explicitly to the distribution is the pnorm(y) in the ad.statistic function. This would be replaced by p___(y), where ___ refers to the relevant CDF function for you, and the parameters would have to be redefined appropriately.