How to generate Random Number for User defined model

Hello everyone,
I have a probability model and I would like to perform a posterior predictive check in Stan using rstan package in R. But I could not define the model in generated quantities block, please help me how to generate random numbers. I have shared my model below
f(x) = \frac{{\lambda \beta }} {{x^2 (1 - e^{ - \lambda } )}}\left\{ {1 + \alpha \left( {1 - e^{ - \beta /x} } \right)} \right\}\exp \left( { - \frac{\beta }{x} - \alpha e^{ - \beta /x} } \right)\exp \left\{ { - \lambda \left( {1 - e^{ - \beta /x} } \right)\exp ( - \alpha e^{ - \beta /x} )} \right\}{\text{ }}

functions{real poissonISG_lpdf(real y, real alpha, real beta, real lambda){
  return log(beta)+log(lambda)- log(1 - exp(-lambda))- 2*log(y) + log(1+alpha*(1 - exp(-(beta*1/y))))- 
  (beta*1/y) - alpha*exp(-(beta*1/y)) -lambda*(1-exp(-(beta*1/y)))*exp(-alpha*exp(-(beta*1/y)));
 }
}

data{
  int N;
  real y[N];
}

parameters{
  real <lower=0> alpha;
  real <lower=0> beta;
  real <lower=0> lambda;
}

model{
  for(i in 1 : N){
    y[i]~ poissonISG(alpha, beta, lambda);
  }
  alpha~ gamma(0.5, 0.5);
  beta~ gamma(0.1, 0.1);
  lambda~ gamma(0.1, 0.1);
}

generated quantities{ please help}

my data set is

y <- c(1.312, 1.314, 1.479, 1.552, 1.700, 1.803, 1.861, 1.865, 1.944, 1.958, 
       1.966, 1.997, 2.006, 2.021, 2.027, 2.055,2.063, 2.098, 2.14, 2.179, 2.224, 2.240, 
       2.253, 2.270, 2.272, 2.274, 2.301, 2.301, 2.359, 2.382, 2.382, 2.426,2.434, 2.435, 
       2.478, 2.490, 2.511, 2.514, 2.535, 2.554, 2.566, 2.57, 2.586, 2.629, 2.633, 2.642, 
       2.648, 2.684,2.697, 2.726, 2.770, 2.773, 2.800, 2.809, 2.818, 2.821, 2.848, 2.88, 
       2.954, 3.012, 3.067, 3.084, 3.090, 3.096,3.128, 3.233, 3.433, 3.585, 3.585)

usually, people do via transformation of the uniform number generator. You may check, first, how to generate the exponential distribution using uniform rng. In more complicated cases as yours, it becomes more like a mathematical problem

ps: see an example RNG for truncated distributions - #6 by bgoodri

Thanks aakhmetz for providing a helpful link, but I could not find the inverse of my CDF as in the link provided by you, example below

real weibull_lb_rng(real alpha, real sigma, real lb) {
    real p = weibull_cdf(lb, alpha, sigma);  // cdf for bounds
    real u = uniform_rng(p, 1);
    return sigma * (-log1m(u))^(1 / alpha);  // inverse cdf for value
  }

I have shared my CDF below
F(x){\text{ }} = {\text{ }}1 - \frac{1} {{(1 - e^{ - \lambda } )}}\left[ {1 - \exp \left\{ { - \lambda \left( {1 - e^{ - \beta /x} } \right)\exp ( - \alpha e^{ - \beta /x} )} \right\}} \right]{\text{ }};\;\alpha \beta \lambda > 0,x > 0
Please help me how to fix these issues.
Thanks

Dear @Laxmi_Pd_Sapkota, as I wrote earlier, it is definitely a mathematical problem, and how you solve it depends on your mathematical skills (personally, I don’t know how to find the analytical form for the inverse CDF of your function)

Is there any function in Stan to obtain the inverse function of CDF like uniroot function in R? Can we draw random numbers using PDF in Stan rather than inverse CDF? Plz share any links or documents if any to fix these issues.

there is algebraic solver, you may check the Stan guide

In the example of algebra solver chapter provides to calculate roots numerically but
my query is could we find a function as a root that can be used to RNG function