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?
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.
-
1$\begingroup$ +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 tooptim
and (2) you can just as easily handle other families of distributions. $\endgroup$ – whuber♦ Aug 1 '12 at 21:02 -
1$\begingroup$ @whuber - No, I didn't realize that my answer was a special case of your suggestion, although I could see a similarity. Now that I've actually thought about it, though, nice! (+1) Another tool in the toolbox. As for the OP's objective, the difficulty I saw is that if he assumes a distributional form, there is information about the tails in the non-tail data, so somehow I feel reluctant to throw the non-tail data out completely. But if it's only a description of the tails of the data, and the distributional assumption is only for convenience, I can see the difference in objectives. $\endgroup$ – jbowman Aug 2 '12 at 14:36
fitdistr
. $\endgroup$ – user10525 Aug 1 '12 at 18:51