randp.m

This function generates Pareto random variables (of type I). See
"Statistical Distributions", Evans, Hastings and Peacock, Wiley, 1993
or http://www.maths.adelaide.edu.au/matthew.roughan/probability_distrns/node6.html
or http://en.wikipedia.org/wiki/Pareto_distribution

The Pareto distribution is a classic "heavy-tailed" or "power-law" distribution. It has distribution function
F(x) = 1 - (b/x)^alpha, for x>=b
and density
f(x) = (alpha/b) * (b/x)^(alpha+1), for x>=b

Its mean is
E[X] = b * alpha/(alpha-1), for alpha>1
but note that the mean is infinite for alpha<=1

Its variance is
Var(X) = b^2 * alpha/[(alpha-1)^2*(alpha-2)], for alpha>1
but note that the variance is infinite for alpha<=2

Median
Med(X) = b * 2^(1/alpha)

Note the function calls "rand" so if you want to control the seed, use rand('state', seed).

Also note that the type II Pareto is just shifted so that x>=0, so to obtain this distributio just take
x = randp(N,M,alpha,b) - b;