How do I fit a set of data to a Pareto distribution in R?

Have, let's say, the following data:

[1] 8232302 684531 116857 89724 82267 75988 63871 23718 1696 436 439 >[12] 248 235

Want a simple way to fit this (and several other datasets) to a Pareto distribution. Ideally it would output the matching theoretical values, less ideally the parameters.

-

Well, if you have a sample $X_1, ..., X_n$ from a pareto distribution with parameters $m>0$ and $\alpha>0$ (where $m$ is the lower bound parameter and $\alpha$ is the shape parameter) the log-likelihood of that sample is:

$$n \log(\alpha) + n \alpha \log(m) - (\alpha+1) \sum_{i=1}^{n} \log(X_i)$$

this is a monotonically increasing in $m$, so the maximizer is the largest value that is consistent with the observed data. Since the parameter $m$ defines the lower bound of the support for the Pareto distribution, the optimum is

$$\hat{m} = \min_{i} X_i$$

which does not depend on $\alpha$. Next, using ordinary calculus tricks, the MLE for $\alpha$ must satisfy

$$\frac{n}{\alpha} + n \log( \hat{m} ) - \sum_{i=1}^{n} \log(X_i) = 0$$

some simple algebra tells us the MLE of $\alpha$ is

$$\hat{\alpha} = \frac{n}{\sum_{i=1}^{n} \log(X_i/\hat{m})}$$

In many important senses (e.g. optimal asymptotic efficiency in that it achieves the Cramer-Rao lower bound), this is the best way to fit data to a Pareto distribution. The R code below calculates the MLE for a given data set,X.

pareto.MLE <- function(X)
{
n <- length(X)
m <- min(X)
a <- n/sum(log(X)-log(m))
return( c(m,a) )
}

# example.
library(VGAM)
set.seed(1)
z = rpareto(1000, 1, 5)
pareto.MLE(z)
[1] 1.000014 5.065213


Edit: Based on the commentary by @cardinal and I below, we can also note that $\hat{\alpha}$ is the reciprocal of the sample mean of the $\log(X_i /\hat{m})$'s, which happen to have an exponential distribution. Therefore, if we have access to software that can fit an exponential distribution (which is more likely, since it seems to arise in many statistical problems), then fitting a Pareto distribution can be accomplished by transforming the data set in this way and fitting it to an exponential distribution on the transformed scale.

-
(+1) We can write things a bit more suggestively by noting that $Y_i = \log(X_i/m)$ is distributed exponential with rate $\alpha$. from this and the invariance of MLEs under transformation we conclude at once that $\hat\alpha = 1/\bar Y$, where we replace $m$ by $\hat m$ in the latter expression. This also hints at how we might use standard software to fit a Pareto even if no explicit option is available. –  cardinal May 1 '12 at 1:52
@cardinal - So, $\hat{\alpha}$ is the reciprocal of the sample mean of the $\log(X_i/\hat{m})$'s, which happen to have an exponential distribution. How does this help us? –  Macro May 1 '12 at 2:27
Hi, Macro. The point I was trying to make was that the problem of estimating the parameters of a Pareto can be (essentially) reduced to that of the estimation of the rate of an exponential: Via the transformation above, we can convert our data and problem into a (perhaps) more familiar one and immediately extract the answer (assuming we, or our software, already know what to do with a sample of exponentials). –  cardinal May 1 '12 at 11:59
How can i measure the error of this kind of fit? –  emanuele Jun 18 '12 at 8:35
@emanuele, the approximate variance of an MLE is the inverse of the fisher information matrix, which will require you to calculate at least one derivative of the log-likelihood. Or, you could use a kind of bootstrap resampling to estimate the standard error. –  Macro Jun 18 '12 at 13:12