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?
 A: 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))

A: 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$.
A: 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)
 curve(1.4*dexp(x), add=T, col="blue")
 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")
 curve(1.4*exp(-abs(x)),  col="maroon", add=T)
 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.
A: 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$.
