How to generate a random number from a pareto distribution I'm working on a problem where I am trying to generate a random number from a Pareto distribution.


*

*Using some measured data,  I have been able to fit a Pareto distribution to this data set with shape/scale values of $4/6820$ using the R library fitdistrplus.

*Now I want to, using the above scale and shape values to generate random numbers from this distribution.

*Doing some searching here (http://www.ntrand.com/pareto-distribution/)  and here (http://www.randomservices.org/random/special/Pareto.html), it mentions that the inverse function or quantile function can be used to generate a random number i.e. 
$$x=\frac{b}{(1−U)^{1/a}} \in U(0,1)$$
Using two methods (some C code and Excel formula) I've generated around $1000$ numbers using my shape & scale parameters. However when I import the data back into R and formally validate whether the randomly generated numbers are from a Pareto distribution an Anderson Darling GoF tests tells me that the Pareto distribution is not a good fit, additionally the fitdist library tells me my shape/scale parameters are totally different than the original parameters. (e.g. shape $= 8260157$, scale $= 10834967$)
Any ideas what I am doing wrong here?
Regards
Jonathan
 A: Either something is wrong with your numerical generation or with your maximum likelihood analysis
If the shape is $a$ and the scale is $b$, I would agree with you that given a uniform random variable $U$ on $(0,1)$ you would have $$X=\dfrac{b}{(1-U)^{1/a}}$$ as having a Pareto distribution.
Given a set of values drawn from a Pareto distribution, I would have thought that the maximum likelihood estimate of the scale would be $$\hat{b} = \min(x_i)$$ and of the shape would be $$\hat{a} = \dfrac{1}{\overline{\log(x_i)} - \log\left(\hat{b}\right)}$$
The following R code suggests this seems to work, since 
shape <- 4
scale <- 6820 
cases <- 1000
set.seed(2016)
u <- runif(cases)
x <- scale / (1-u)^(1/shape) 
mlescale <- min(x)
mleshape <- 1 /( mean(log(x)) - log(mlescale) )

gives
> mleshape 
[1] 3.968604
> mlescale 
[1] 6824.48

which are close to to original values 
A: If $$U\in [0,1)$$ is uniformly distributed random variable, then $$E(\lambda)=-\ln(1-U)/\lambda$$ is exponentially distributed random variable, and then $$P(x_m,\alpha)=x_m\cdot\exp(E(\alpha))$$ is Pareto distributed random variable with scale $x_m$ and shape $\alpha$. $$L(x_m,\alpha)=P(x_m,\alpha)-x_m$$ generates Lomax distributed random variable.
A: Ok i believe i figured out the answer, it looks like the formulae that i mentioned in my original post was not correct. I figured out that the numbers being generated where out by a factor the the scale i used i.e. 6820.53374 
Apologies I am using my c code to show the correct computations (See below):
double pareto_arrival(void) {
        /* R computations estimate that Pareto is a reasonable fit with /
        / shape = 4.14104, scale = 6820.53374 paramaters /
        / https://cran.r-project.org/web/packages/actuar/actuar.pdf */
float a;
float b;
float inv_fun_denom;
float par_arr;

a =  4.14104;
b = 6820.53374;

inv_fun_denom = pow(1-drand48(), 1/a);
par_arr = (b/inv_fun_denom)-b; //adding the -b did the trick

return par_arr;

}
Page 64 of this pdf helped put me on the right track 
https://cran.r-project.org/web/packages/actuar/actuar.pdf
