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.


2 Answers 2


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.

  • 3
    $\begingroup$ (+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. $\endgroup$
    – cardinal
    May 1, 2012 at 1:52
  • $\begingroup$ @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? $\endgroup$
    – Macro
    May 1, 2012 at 2:27
  • 2
    $\begingroup$ 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). $\endgroup$
    – cardinal
    May 1, 2012 at 11:59
  • $\begingroup$ How can i measure the error of this kind of fit? $\endgroup$
    – emanuele
    Jun 18, 2012 at 8:35
  • $\begingroup$ @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. $\endgroup$
    – Macro
    Jun 18, 2012 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
  • $\begingroup$ Should that be library(fitdistrplus) ? $\endgroup$
    – Sean
    Jan 27, 2017 at 13:57
  • 1
    $\begingroup$ @Sean yes, editing response accordingly $\endgroup$ May 9, 2017 at 1:07
  • $\begingroup$ Note that the call to library(actuar) is required for this to work. $\endgroup$
    – jsta
    Dec 14, 2017 at 21:08
  • $\begingroup$ What does fp$estimate["shape"] represent in this case? Is it perhaps the estimated alpha? Or beta? $\endgroup$ Dec 17, 2018 at 20:39
  • $\begingroup$ actuar package in not available in R $\endgroup$
    – 89_Simple
    Jul 19, 2021 at 23:24

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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