# 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);
}



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