# Fitting the tails of the distribution

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?

-
"best" in what sense? Note also that any method based on fitting the tails would require a moderate-large sample size. If you just want to fit the distribution using MLE, then you can use the command fitdistr. – user10525 Aug 1 '12 at 18:51
@Procrastinator - I am assuming that a normal distribution is the "best fit" distribution for the data. The data consists of 300 data points. I used the command 'fitdistr' using MLE, but I believe this particular function matches the entire data. I only want to focus in on the tails. – Hank Aug 1 '12 at 18:56
How are you defining the tails? Are you interested on estimating "optimally" ${\mathbb P}(\vert X \vert > T)$, with $X\sim N(\mu,\sigma)$, for a certain $T$? – user10525 Aug 1 '12 at 18:59
A general method, somewhat ad hoc but very effective, is to use weighted regression of the probability plot, focusing the weights in the tails. – whuber Aug 1 '12 at 19:17
There is a large literature on extreme value theory that attempts to model the tails of a distribution. This package has some tools that can be of help cran.r-project.org/web/packages/evir/index.html – John Aug 1 '12 at 19:27

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)
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)
[1] 0.6334132
... 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")

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.
+1 It occurs to me that your sense of "best" fit might not be the same as the OP's, who might only be seeking a parsomonious description of the data that gives a reliable picture of the tails. In that case, there are no small-sample concerns and asymptotic convergence is irrelevant (although nice). BTW, did you notice that your answer is a special case of the probability plotting approach I suggested? The import of that is (1) you can easily compute the AD statistic using, say, lm without recourse to optim and (2) you can just as easily handle other families of distributions. – whuber Aug 1 '12 at 21:02