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?

  • 1
    $\begingroup$ "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. $\endgroup$ – user10525 Aug 1 '12 at 18:51
  • $\begingroup$ @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. $\endgroup$ – Hank Aug 1 '12 at 18:56
  • 1
    $\begingroup$ 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$? $\endgroup$ – user10525 Aug 1 '12 at 18:59
  • 2
    $\begingroup$ A general method, somewhat ad hoc but very effective, is to use weighted regression of the probability plot, focusing the weights in the tails. $\endgroup$ – whuber Aug 1 '12 at 19:17
  • 1
    $\begingroup$ @vdesai The probability plot approach works for scale-location families like the Normal because changes of scale and location correspond to changes in slope and intercept (respectively) in the probability plot. $\endgroup$ – whuber Aug 2 '19 at 12:50

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")

enter image description here

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(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")
AD:   -0.1817909    1.145518 
MLE:  -0.2015397    1.091143 

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

enter image description here

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 to optim 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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.