Finding random number based on a specific discrete distribution using R I have a discrete distribution with cumulative function:
$$
F(x)=Pr(X \leq x) = \frac{1}{(x-1)!}
\left( \Gamma(x,\lambda)-\frac{\lambda^{x-1} \exp(-\lambda)}{\lambda+1} \right)
$$
for $x=1,2,3,\ldots . $
If I am not wrong, if $U\sim \textrm{unif}(0,1)$, then I can get the random data using $x=F^{-1}(u)$. But it is not easy to get the inverse of the above $F(x)$. My goal is to generate data which follows a distribution with cumulative function $F$. How do I continue with it? I tried the following:
pusbp = function(x){
  f1 = factorial(x-1)
  f2 = gammainc(x,m)
  f3 = (m^x)*exp(-m)/(m+1)
  
  Fx = (f2-f3)/f1
  Fx
}

inverse.cdf<-function(x,pusbp,starting.value=0){
  lower.found<-FALSE
  lower<-starting.value
  while(!lower.found){
    if(pusbp(lower)>=(x-.000001))
      lower<-lower-(lower-starting.value)^2-1
    else
      lower.found<-TRUE
  }
  upper.found<-FALSE
  upper<-starting.value
  while(!upper.found){
    if(pusbp(upper)<=(x+.000001))
      upper<-upper+(upper-starting.value)^2+1
    else
      upper.found<-TRUE
  }
  uniroot(function(y) pusbp(y)-x,c(lower,upper))$root
}


vars<-apply(matrix(runif(1000)),1,function(x) inverse.cdf(x,pusbp))
hist(vars)

This is a copy paste code with slight tweak from here. I don't understand the codes.
The error I get is:
Error in if (pusbp(lower) >= (x - 1e-06)) lower <- lower - (lower - starting.value)^2 -  : 
  missing value where TRUE/FALSE needed
In addition: Warning message:
In gamma(x + 1) :
 Error in if (pusbp(lower) >= (x - 1e-06)) lower <- lower - (lower - starting.value)^2 - : 
missing value where TRUE/FALSE needed


Alternate approach to computing the inverse cdf is to use the fact that we can generate a draw as $U\sim Unif(0,1)$ then $X=\underset{i}{\textrm{min}} \{i :  P(X\leq i)\geq U\}$.
Here's an example:
N_max <- 10
x <- seq(1,N_max)
m <- 1
fx <- sapply(x, pusbp)

We start by sampling from the distribution with $\texttt{m}=1$ with $\texttt{m}$ denoting $\lambda$. At this value of $\lambda$, $P(X\leq 10) = 0.9999994$, so we aren't missing much by truncating this RV at 10. More on that later. But here's a basic implementation of a sampler for that distribution:
n_samp <- 1000
u <- runif(n_samp)
comp <- outer(u,fx,function(x,y) x<y)
samp <- apply(comp, 1, function(x) min(which(x)))
  

Which gives us an empirical CDF close to nominal (open circles denote theoretical values):
plot(ecdf(samp))
points(x,fx)


To make this rigorous let's circle back to choosing $\texttt{N_max}$, the number of points to evaluate the cdf at. If we want to make sure that we never truncate the variable to be smaller than it should be, we can draw $u_j \overset{iid}{\sim} U(0,1)$ first, and then make sure to calculate $F(x)$ out to at least $\underset{i}{\textrm{min}} \{i :  P(X\leq i)\geq \max u_j \}$.
