# how to generate data from cdf which is not in closed form?

i am working on a distribution whose pdf and cdf is $$f(x,\alpha,\beta)=\frac{(\frac{\beta}{\alpha})(\frac{x}{\alpha})^{\beta}}{(1+(\frac{x}{\alpha})^{\beta})^{2}}\frac{\sin(\frac{\pi}{\beta})}{(\frac{\pi}{\beta})}$$ $$x>0$$, $$\alpha>0$$, $$\beta>0$$ $$F(x)= \frac{\sin(\frac{\pi}{\beta})}{(\frac{\pi}{\beta})}\frac{1}{\beta}(\frac{x}{\alpha})^{1-\beta}\log{(1+(\frac{x}{\alpha})^{\beta})}-\frac{(\frac{x}{\alpha})}{1+(\frac{x}{\alpha})^{\beta})}-(\frac{1-\beta}{\beta})[(\frac{x}{\alpha})+\sum_{r=1}^{\infty}\frac{(-1)^{r}(\frac{x}{\alpha})^{1+r\beta}}{r(1+r\beta)}]$$ clearly, cdf is not in neat form. So how do I generate data from this cdf in r as inverse cdf method wont't apply here?

• Please provide complete definitions. Your $F$ is not a CDF; evidently its argument $x$ is intended to lie within some interval (it looks like $[0,1]$ but I don't want to have to verify that). What are the possible values of $(\alpha,\beta)$? – whuber Feb 5 at 18:14
• There's no guarantee $F$ is a CDF unless it is correctly normalized. Re the edits: from $f$ it looks like $x^\beta$ has an F distribution, but your expression for its CDF is a strange one. It won't even converge for $x\gt \alpha.$ Where does it come from? Are you sure there aren't any typographical errors in your representation of $f$? – whuber Feb 5 at 18:21
• There must be typos. For instance, $f$ is directly proportional to $\alpha^{1-\beta},$ which means it can be a valid PDF for at most one value of $\alpha$--but then $\alpha$ would not be a parameter. – whuber Feb 5 at 18:43
• It still doesn't work, for the same reason. Most likely the "$x/\beta$" in the denominator should be "$x/\alpha.$" However, it's important that you get this right. There are still obvious typos. For instance, since $\sin(\pi/\beta)$ can be negative, it cannot possibly be a factor either of $f$ or $F.$ One way would be to describe the distribution in detail so people could derive it themselves. Another would be to quote an accessible reference. – whuber Feb 5 at 19:06
• Reposting your question did not fix the fact that you haven't described a valid distribution. Please visit our help center for information on how this site works. – whuber Feb 6 at 21:12

### Solution

By exploiting the principles explained at https://stats.stackexchange.com/a/495347/919, we can examine the expression for

$$f_X(x,\beta,\alpha) \propto \frac{(x/\alpha)^{\beta+1}}{\left(1 + (x/\alpha)^\beta\right)^2}\frac{\mathrm{d}x}{x}$$

and see that if it is to make any sense at all (a question I will prescind from until later), then any random variable $$X$$ with this density must be the $$1/\beta$$ power of some random variable $$Z$$ with an F ratio distribution. The latter's density is determined by two parameters (its "degrees of freedom," both of which are positive) as

$$f_Z(z,\nu_1,\nu_2) \propto \frac{z^{\nu_1/2}}{\left(1 + \frac{\nu_1}{\nu_2}z\right)^{(\nu_1+\nu_2)/2}} \frac{\mathrm{d}z}{z}.$$

Setting $$(x/\alpha)^\beta = z$$ and comparing these expressions gives

\left\{\begin{aligned}\nu_1 &= 2 + \frac{2}{\beta},\\\nu_2 &= 2 - \frac{2}{\beta}.\end{aligned}\right.

Since $$\nu_2\gt 0,$$ this works only if $$\beta \gt 1.$$ At the end of this post I will demonstrate that $$f_X$$ is undefined for $$0\lt \beta\le 1,$$ showing that nothing is lost by viewing $$X$$ as a power of an $$F_{2+2/\beta,2-2/\beta}$$ variable.

(Incidentally, the series expansion in the question fails to converge for $$x\gt \alpha$$ and so is invalid. The CDF of an F-ratio distribution is given by a regularized incomplete Beta function.)

Consequently, (1) $$f_X$$ does define a distribution for $$\beta \gt 1$$ and (2) to draw a random variate $$X$$ with this distribution,

Draw a random variable $$Z$$ from an F-ratio distribution with $$2+2/\beta$$ and $$2-2\beta$$ degrees of freedom and set $$X = \alpha \left(Z \left(\frac{\beta+1}{\beta-1}\right)\right)^{1/\beta}.$$

### Practical matters

Because the F-ratio distribution is widely used in statistics and is simply related to other common distributions (such as Beta, Gamma, and Chi-squared distributions), efficient, accurate methods to compute with it abound. To generate random variates with $$F$$ distribution I used the rf function in R.

As a demonstration, for three widely varying values of $$\beta$$ I generated 100,000 independent realizations of $$X$$ in this manner, plotted the histogram of $$\log(X)$$ (because the distribution can be very skewed), and overplotted that with the graph of $$f_X(;\beta,\alpha).$$ The agreements are excellent.

(The vertical lines shows the modes of these distributions.) The R command to produce these data was

 x <- alpha * (rf(n, 2+2/beta, 2-2/beta) * (beta+1)/(beta-1))^(1/beta)


### Further analysis and limitations

Let's derive the distribution given by $$f_X.$$ Suppose $$Y$$ has a standard logistic distribution, which means it has a density $$F_Y$$ defined on all real numbers $$y$$ by

$$f_Y(y) = \frac{e^y}{(1+e^y)^2} = \frac{\mathrm{d}}{\mathrm{d}y} \frac{1}{1 + e^{-y}}.$$

Because the right hand side (the antiderivative) increases from $$0$$ to $$1,$$ it is a cumulative distribution function, whence $$f_Y$$ is a probability density function.

Now suppose $$Y$$ is the logarithm of some (positive) random variable $$X;$$ that is, $$X = e^Y.$$ Then the probability element of $$X$$ is

$$f_X(x)\mathrm{d}x = f_Y(\log x)\mathrm{d}(\log x) = \frac{x}{(1+x)^2}\frac{\mathrm{d}x}{x} = \frac{1}{(1+x)^2}\mathrm{d}x.$$

When $$Y$$ is rescaled by a factor $$\beta\gt 0,$$ the rules of exponentiation imply $$X$$ is transformed to $$X^\beta.$$ Under this transformation the probability element $$\mathrm{d}x/x$$ is merely multiplied by $$\beta,$$ showing immediately that the density of $$X^\beta$$ is proportional to

$$f_X(x,\beta) \propto \frac{x^\beta}{\left(1 + x^\beta\right)^2} \frac{\mathrm{d}(x^\beta)}{x^\beta} \propto \frac{x^{\beta-1}}{\left(1 + x^\beta\right)^2}\mathrm{d}x.$$

Introducing a scale factor $$\alpha$$ for $$X$$ (which would be the exponential of a location parameter for $$Y$$) gives the log-logistic distribution family,

$$f_X(x,\beta,\alpha) \propto \frac{(x/\alpha)^{\beta-1}}{\left(1 + (x/\alpha)^\beta\right)^2}\mathrm{d}x.$$

This distribution can be "length extended" in various ways. By comparison to the question, it is evident that $$f_X(x,\beta,\alpha)$$ has been changed to $$xf_X(x,\beta,\alpha).$$ This weights the probability densities at larger $$x\gt 0$$ directly proportional to $$x,$$ skewing the distribution to the right. This helps us understand what $$f_X$$ is, what it is intended to accomplish, and how it relates to familiar distributions.

It is time to consider the implicit proportionality constants. For this length-extended log-logistic (LELL) distribution to be defined, it must be possible to scale all densities so they integrate to unity. The integral of $$x f_X(x,\beta,\alpha)$$ is (by definition) the expectation of $$X,$$ equal to

$$E[X] = \frac{\alpha \pi/\beta}{\sin(\pi/\beta)}$$

provided $$\beta\gt 1.$$ This is why its reciprocal appears as a factor of $$f_X$$ in the question.

When $$0 \lt \beta \le 1,$$ the LELL construction does not work. Let's see why not. The scale factor doesn't matter, so take $$\alpha=1.$$ For any $$T\gt 1,$$ a series of elementary estimates based on $$x^\beta \gt 1$$ and $$x^{-\beta}\ge x^{-1}$$ when $$x\gt 1$$ (when $$0\lt \beta\le 1$$) gives

\begin{aligned} E[X\mid \beta,1] &= \int_\mathbb{R} x f_X(x,\beta,1)\,\mathrm{d}x\\ &= \int_0^\infty x \frac{x^{\beta-1}}{\left(1 + x^\beta\right)^2}\mathrm{d}x \\ &\ge \int_1^\infty \frac{x^{\beta}}{\left(1 + x^\beta\right)^2}\mathrm{d}x \\ &\gt \int_1^\infty \frac{x^{\beta}}{\left(x^\beta + x^\beta\right)^2}\mathrm{d}x \\ &= \frac{1}{4}\int_1^\infty x^{-\beta}\,\mathrm{d}x\\ &\gt \frac{1}{4}\int_1^T x^{-\beta}\,\mathrm{d}x\\ &= \frac{1}{4(1-\beta)}\left(T^{1-\beta}-1\right). \end{aligned}

Because the final expression grows arbitrarily large as $$T$$ grows, $$E[X]$$ cannot be finite when $$\beta \le 1.$$ This makes it impossible to construct an LELL distribution for $$\beta \le 1.$$

### Working code

This is a series of R commands. On my system, the command to draw independent values following the $$f_X$$ density will produce four million (4,000,000) values per second (on one core).

#
# The PDF in the question.
#
f <- function(x, beta, alpha=1) {
y <- (x/alpha)^beta
beta^2 * sin(pi/beta) / pi * y / (1 + y)^2 / alpha
}
#
# The mode of log(X) where the density of X is f.
#
mode <- function(beta, alpha=1) {
log(alpha) + log((1+beta)/(beta-1)) / beta
}
#
# Simulations.
#
n <- 1e5                # Simulation size
log.alpha <- 1
alpha <- exp(log.alpha) # Scale; must be positive
beta <- c(1.1, 2, 1e3)  # The beta values to use

par(mfrow=c(1,length(beta)))
set.seed(17)
for (beta in beta) {
#
# Generate random variates.
#
x <- alpha * (rf(n, 2+2/beta, 2-2/beta) * (beta+1)/(beta-1))^(1/beta)
#
# Plot a histogram of log(x) (because X is so skewed).
#
m <- mode(beta, alpha)
ymax <- f(exp(m), beta, alpha)*exp(m)
hist(log(x), freq=FALSE, breaks=80, ylim=c(0, ymax), col="#f0f0f0", border="Gray")
abline(v = m) # Mark the mode
mtext(bquote(list(alpha==e^.(log.alpha), beta==.(beta))), side=3, line=0, cex=0.7)
#
# Overplot the PDF to check.
#
curve(f(exp(x), beta, alpha)*exp(x), add=TRUE, n=501, lwd=2, col="Red")
}
par(mfrow=c(1,1))

• Thank you for this nice explanation. Please clear one more thing to me. I have to generate data from this distribution using Adaptive progressive Type II censoring. In last step of data generation under this censoring, I have to generate data from truncated distribution $$\frac{f(x)}{1-F(x_{j})}$$ where $F(x_{j})$ is value of CDF at $x_{j}$. So is it correct to generate \textbf{Z} from truncated F distribution using rtrunc function in R and then transform it into $X$ using the above method that you explained. – Pulkit Srivastava Feb 9 at 20:49
• I can't answer that because R does not include an rtrunc function. Since you now have an account of the distribution function and it can be inverted (use qf in R), to generate samples from the distribution truncated to the interval $[x_j,\infty)$ you can apply the method I have just described at stats.stackexchange.com/a/508780/919 by taking $a=x_j$ and $b=\infty.$ – whuber Feb 9 at 20:53
• rtrunc function is used with the help of trucndist package in R . Here my f(x) is same LELL distribution whose data generation you explained above. – Pulkit Srivastava Feb 10 at 5:29
• So far I understand, I should generate samples from truncated F ratio distribution and then convert it into X using $$X= \alpha(Z(\frac{\beta+1}{\beta-1}))^\frac{1}{\beta}$$ where Z is from truncated F ratio distribution – Pulkit Srivastava Feb 10 at 5:35

### Searching for something related

• Step 1 use a transformation to simplify te distribution. With $$t = (x/\alpha)^\beta$$ or $$x=\alpha t^{1/\beta}$$ it becomes

$$f(t) \propto \frac{t^{1/\beta}}{(1+t)^{2}}$$

• Step 2. Compare with other distributions. We can use the following quotient to see if this is a Pearson distribution

$$f'(t)/f(t) = \frac{\frac{1}{1-2 \beta}+t}{(t+t^2)/(1/\beta-2)}$$

which is indeed like a Pearson distribution.

• Searching a bit further we see that it is related to the type VI distribution, or F-distribution of Beta prime distribution. The Beta prime distribution below resembles directly with what we had above

$$\frac{1}{B(a,b)} \frac{x^{a-1}}{(1+x)^{a+b}}$$

this works if your $$1/\beta<1$$ (because the beta prime distribution has the restriction $$a>0$$ and $$b>0$$)

So you can sample $$t$$ from the beta prime distribution with any of the software packages for this function, and then transform to $$x$$.

I cannot make sense of your PDF. So I will show an acceptance-rejection method using a PDF that does make sense to me.

The method is similar to that of @Ertxiem (+1), except that I guess mine works better if your (correct) density function has support $$(0,\infty).$$ Both methods are general ones, and I hope you can use one of them to sample from your density when you get its proper formulation.

Suppose you want to sample from the 'half normal' distribution of $$X =|Z|,$$ where $$Z$$ is standard normal, so that the density function of $$X$$ is $$f(x) = \frac{2}{\sqrt{2\pi}}e^{0.5x^2},$$ for $$x > 0.$$ One can show that $$E(X) = \sqrt{2/\pi},\;$$ $$Var(X) = 1 - 2/\pi.$$ We want to take a sample of size $$n = 10^6$$ from this distribution. We cannot write the CDF of this density function or its inverse in closed form.

For this demonstration, forget that the standard normal CDF has been extensively tabled and that its inverse function is implemented in R and other statistical software.

One solution is to find an 'envelope function' $$g(x)\ge f(x),$$ for $$x > 0,$$ where $$g(x)$$ is a multiple of the density function of a distribution from which we do know how to sample. Here, we can take $$g(x) = 1.4e^{-x},$$ which is $$1.4$$ times the density function of $$\mathsf{Exp}(1).$$

hdr = "Half Normal Density with Envelope"
curve(2*dnorm(x), 0, 4, ylim=c(0,1.4), lwd=2, ylab="Density", main=hdr)
abline(h=0, col="green2"); abline(v=0, col="green2")


Now we generate a random sample y of size $$n$$ from $$\mathsf{Exp}(1),$$ using its quantile function and standard normal random variables.

n = 10^6;  y = -log(runif(n)) # exponential
u = runif(n, 0, 1.4*exp(-y))
acc = u < 2*dnorm(y)
x = y[acc]

# for verification: about 2 or 3 places
mean(x);  sqrt(2/pi)
[1] 0.7974243  # aprx E(X)
[1] 0.7978846  # exact
sd(x);  sqrt(1 - 2/pi)
[1] 0.6024202   # aprx SD(X)
[1] 0.6028103   # exact

hdr = "Historam of Sampled Points with Density"
hist(x, prob=T, col="skyblue2", main=hdr)
curve(2*dnorm(x), col="maroon", 0, 5, add=T)


The following figure shows some of the accepted points (blue) and some of those not accepted (orange). [For clarity in the figure, not all simulated points are shown.] The vector x contain the x-coordinates of all accepted points.

R code for image above:

curve(2*dnorm(x), 0, 4, lwd=2, col="blue", ylim=c(0,1.42), xaxs="i",
ylab="Density", main="Accepted (x) Beneath Folded Normal Density")
abline(h=0, col="green2")
plt = 20000; y.pl = y[1:plt] # plot fewer points for clarity
u.pl = u[1:plt];  acc.pl=acc[1:plt]
points(y.pl, u.pl, pch=".", col="orange")
points(y.pl[acc.pl], u.pl[acc.pl], pch=".", col="skyblue4")


Notes: (1) In the first block of code to generate values x, I could have used rexp(n) to generate y instead of using the quantile function of uniform random variables with -log(runif(n)). Also, instead of 2*dnorm(y), I could have used (2/sqrt(2*pi))*exp(-0.5*y^2).

(2) I'm not sure if you'll know the mean, variance, or other facts about your density function, so you may not have anything to 'verify' after taking your sample x. But you could still plot a histogram of the simulated data and overlay your density function.

(3) In order for my method to work, you will have to be able to find some distribution that can be simulated, and for which a multiple of the density majorizes your density.

• I'm pretty sure $\beta \gt 1.$ In that case, assuming the expression for $f$ in the question is correct (and it would then make sense), $(x/\alpha)^\beta$ should have an $F$ distribution and the problem is solved. For reasons I gave in comments to the question, it might be best to ignore the expression for the CDF. – whuber Feb 6 at 23:34
• @whuber. If your comment refers to my Answer, I can't find the 'expression for the CDF' to which you refer. Glad to fix anything misleading. – BruceET Feb 7 at 0:37
• @whuber could you please explain why f is not a pdf. It is pdf of length baised log logistic distribution. Research paper regarding this have been published already. – Pulkit Srivastava Feb 7 at 7:18
• Pulkit, in my answer I have included a section showing what the allowable ranges of the parameters are and explaining why that is. – whuber Feb 7 at 15:39
• @whuber. So, properly formulated and parameterized, the distribution has known PDF, CDF, and quantile functions, and so fancy Monte Carlo sampling methods are not required. Thanks. – BruceET Feb 8 at 7:13

You may try to use a Monte-Carlo method.

Generate a random value for $$(x, \alpha, \beta)$$ and $$f_e$$. If $$f_e < f(x, \alpha, \beta)$$ accept that triple, otherwise, reject it.

This is easier to make work if $$(x, \alpha, \beta)$$ and $$f$$ are bounded than if they're unbounded.

Edit:

I will assume that $$x \in [x_L; x_u]$$, $$\alpha \in [\alpha_L; \alpha_U]$$, $$\beta \in [\beta_L; \beta_U]$$ and $$f \in [f_L; f_U]$$, with $$f_L = 0$$.

Step 1: draw randomly from an uniform distribution $$x_e \in [x_L; x_u]$$, $$\alpha_e \in [\alpha_L; \alpha_U]$$, $$\beta_e \in [\beta_L; \beta_U]$$ and $$f_e \in [f_L; f_U]$$.

Step 2: If $$f_e < f(x_e, \alpha_e, \beta_e)$$, then the values $$(x_e, \alpha_e, \beta_e)$$ are valid and you can add that triple to the generated dataset; otherwise reject the triple.

Repeat these two steps until you have a large enough dataset.

Note: in case $$\alpha$$ and $$\beta$$ are parameters, you'll just have to draw $$x_e$$.

• Thanks for your response. But it will be helpful if you could explain more. – Pulkit Srivastava Feb 6 at 14:04
• I am using different R code for generating samples from this dist but not able to justify each steps clearly alpha=1.5 beta=3.2 pdf=function(x)((beta/alpha)*((x/alpha)^beta)*sin(pi/beta))/(((1+(x/alpha)^beta)^2)*(pi/beta)) rand_smplfunc <- function(func, lower, upper){ x_values <- seq(lower,upper,by = 10^-3) sample(x = x_values, size = 1, replace = TRUE,prob = func(x_values)) } rand_smplfunc(pdf,0,10)) – Pulkit Srivastava Feb 6 at 19:57
• Neither this nor any other method will work, because the $f$ is not a PDF and $F$ doesn't even exist when $x\ge \alpha.$ See the discussion in the duplicate thread. – whuber Feb 6 at 21:13
• @whuber could you please explain why f is not a pdf. It is pdf of length baised log logistic distribution. Research paper regarding this have been published already. – Pulkit Srivastava Feb 7 at 7:18