How to simulate a random variable with this density? I need to generate random variables $R$ with this distribution function:
$$F_R(y)=\frac{1}{\mu(\lambda,\tau)}\int_{0}^{y}e^{-\lambda x^{\tau}}dx$$
where $0\lt\tau\lt 1,$ $\lambda\gt 0$ (it is like the rate of an exponential distribution), and $\mu$ is the normalizing constant.
But, I am not clear how to proceed. Is it ok if I consider the associated density; i.e.,
$$f_R(x)=\frac{1}{\mu}e^{-\lambda x^\tau}$$
How to do that?
 A: Clearly $R \ge 0$ (for otherwise the integrand would not be defined when $y \lt 0$ for many such values of $\tau$).  Because the differential element is proportional to
$$f_R(r)\mathrm{d}r \propto \exp(-\lambda  r^\tau) \mathrm{d}r,$$
the variable $U = \lambda R^\tau$ has differential element proportional to
$$f_U(u)\mathrm{d}u \propto u^{1/\tau}\exp(-u) \frac{\mathrm{d}u}{u},$$
exhibiting $U$ as a $\Gamma\left(\frac{1}{\tau}\right)$ variable.  (Thus, $R$ follows a Generalized Gamma Distribution.)  Consequently, you may generate $U$ according to this Gamma distribution (with many well-known efficient methods of generation) and compute
$$R = \left(\frac{U}{\lambda}\right)^{1/\tau}.$$
To illustrate, here is the histogram of $10^4$ independent realizations of $R$ computed in this manner with the R code below.  It is shown on a logarithmic axis for clarity, on which the density of $\log R$ implied by $f_R$ is plotted as a red curve.  The correspondence is excellent:

This eight-year-old workstation generates five million values of $R$ per second.  You might be able to do better with a more modern one.
#-- Specify parameters.
lambda <- 3
tau <- 1/4

#-- Generate variates.
n <- 1e4
R <- (rgamma(n, 1/tau, 1) / lambda)^(1/tau)

#-- Plot the results and compare to the intended distribution.
hist(log(R), freq=FALSE)
curve(exp(x - lambda * exp(tau * x) - lgamma(1+1/tau) + log(lambda)/tau),
      add=TRUE, col="Red", lwd=2)

