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

Question

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.

Answer 1

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.

Answer 2

You can use the fitdist function provided in fitdistrplus package:

library(MASS)
library(fitdistrplus)
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