# How to know if my data fits Pareto distribution?

I have a sample which is a vector with 220 numbers. Here is a link to a histogram of my data.. And I wish to check if my data fits a Pareto distribution, but I don't want to see QQ plots with that distribution, but I need an exact answer with p-value in R, such as Anderson-Darling test for normality (ad.test). How could I do that? Please be as specific as you can.

• The result of a statistical test will not tell you that your data have a Pareto distribution. In fact you can be pretty certain that if it's real data, they don't have a Pareto distribution. All a test will show you is whether you have enough data to pick up the amount deviation from being Pareto that you have. That is, if it rejects all it says is 'yeah the sample size was big enough to tell you what you already knew'. Why would you undertake such an exercise, one which cannot answer the actual question you have? Nov 30, 2013 at 23:14
• Is your question really nothing more than 'what lines of I code do I write to make program R do procedure X'? Then it's off topic here. It might qualify as a programming question. If there's a statistical aspect to your question (like 'does this make any sense to do?') then you should clarify and emphasize those aspects of it Nov 30, 2013 at 23:14
• Now to the Anderson-Darling test (or, for that matter, the Kolmogorov-Smirnov that @Zen suggested above). Those are tests for completely specified distributions. That is, for the tests to have the desired properties you must specify a priori (NOT estimate) all the parameters. So you can't use either of them for this exercise because you don't have prespecified parameters. (Presumably you're doing this at the suggestion of someone else. It's very hard to explain misconceptions to someone via an intermediary.) Nov 30, 2013 at 23:16
• What are you doing this testing for? e.g. what actions will change depending on whether you reject or fail to reject? Nov 30, 2013 at 23:25
• You should always look at a QQ plot, regardless of your motive. And you should not fetishize an "exact" P-value. A different test would give you a different "exact" P-value. Dec 1, 2013 at 9:34

(PS) First of all I think Glen_b is right in his above comments on the usefulness of such a test: real data are surely not exactly Pareto distributed, and for most practical applications the question would be "how good is the Pareto approximation?" – and the QQ plot is a good way to show the quality of such an approximation.

Any way you can do your test with the Kolmogorov-Smirnov statistic, after estimating the parameters by maximum likelihood. This parameter estimation prevents to use the $p$-value from ks.test, so you can do parametric bootstrap to estimate it. As Glen_b tells in the comment, this can be connected to Lilliefors test.

Here are a few lines of R code.

First define the basic functions to deal with Pareto distributions.

# distribution, cdf, quantile and random functions for Pareto distributions
dpareto <- function(x, xm, alpha) ifelse(x > xm , alpha*xm**alpha/(x**(alpha+1)), 0)
ppareto <- function(q, xm, alpha) ifelse(q > xm , 1 - (xm/q)**alpha, 0 )
qpareto <- function(p, xm, alpha) ifelse(p < 0 | p > 1, NaN, xm*(1-p)**(-1/alpha))
rpareto <- function(n, xm, alpha) qpareto(runif(n), xm, alpha)


The following function computes the MLE of the parameters (justifications in Wikipedia).

pareto.mle <- function(x)
{
xm <- min(x)
alpha <- length(x)/(sum(log(x))-length(x)*log(xm))
return( list(xm = xm, alpha = alpha))
}


And this functions compute the KS statistic, and uses parametric bootstrap to estimate the $p$-value.

pareto.test <- function(x, B = 1e3)
{
a <- pareto.mle(x)

# KS statistic
D <- ks.test(x, function(q) ppareto(q, a$xm, a$alpha))$statistic # estimating p value with parametric bootstrap B <- 1e5 n <- length(x) emp.D <- numeric(B) for(b in 1:B) { xx <- rpareto(n, a$xm, a$alpha); aa <- pareto.mle(xx) emp.D[b] <- ks.test(xx, function(q) ppareto(q, aa$xm, aa$alpha))$statistic
}

return(list(xm = a$xm, alpha = a$alpha, D = D, p = sum(emp.D > D)/B))
}


Now, for example, a sample coming from a Pareto distribution:

> # generating 100 values from Pareto distribution
> x <- rpareto(100, 0.5, 2)
> pareto.test(x)
$xm  0.5007593$alpha
 2.080203

$D D 0.06020594$p
 0.69787


...and from a $\chi^2(2)$:

> # generating 100 values from chi square distribution
> x <- rchisq(100, df=2)
> pareto.test(x)
$xm  0.01015107$alpha
 0.2116619

$D D 0.4002694$p
 0


Note that I do not claim that this test is unbiased: when the sample is small, some bias can exist. The parametric bootstrap doesn’t take well into account the uncertainty on the parameter estimation (think to what would happen when using this strategy to test naively if the mean of some normal variable with unknown variance is zero).

PS Wikipedia says a few words about this. Here are two other questions for which a similar strategy was suggested: Goodness of fit test for a mixture, goodness of fit test for a gamma distribution.

• When you adjust the distribution of the test statistic for the estimation of parameters in this way, it's not a K-S test (even though based on a K-S statistics) - it's a particular type of Lilliefors test. This is no longer nonparametric, but one can construct one via simulation for any given distribution. Lilliefors did this specifically for the normal and exponential ... back in the 1960s. Nov 30, 2013 at 23:38
• Thanks for this comment @Glen_b I didn't know that. Dec 1, 2013 at 0:06
• No problem; it changes nothing about the content of what you're doing (which is fine as it is), only what it ought to be called. Dec 1, 2013 at 0:26
• @Glen_b I made some substantial changes in my answer, thanks again! Dec 1, 2013 at 21:27