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

Have, let's say, the following data:

8232302  684531  116857   89724   82267   75988   63871
23718    1696     436     439     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.

-
What is intended by "matching theoretical values"? The expectations of the order statistics given the parameter estimates? Or something else? – Glen_b Oct 7 '15 at 22:55

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

you can use the fitdist function provided in "fitdistrplus" package

library("MASS")
library("fitdistr")
library("actuar")

#suppose data is in dataPar list
fp <- fitdist(dataPar, "pareto", start=list(shape = 1, scale = 500))
#the mle parameters will be stored in fp\$estimate


hope it helps

-