Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

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.

share|improve this question
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
up vote 20 down vote accepted

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. 
z = rpareto(1000, 1, 5) 
[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.

share|improve this answer
(+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


#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

share|improve this answer

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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