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

 A: 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 \}$.
A: Algorithms for discrete lower-bounded distribution
Firstly, for programming these types of functions it is generally best to work in log-probability space, so you should use the log-CDF instead of the CDF.  To examine the general case of this problem, suppose we have a discrete random variable with a programmed log-CDF $L$ over the support $x = x_0,x_0+1,x_0+2,...$ (where the lower bound $x_0$ is known and there is no upper bound) and we want to program the corresponding quantile function and random-generation function for the distribution.
Below I give general algorithms that will compute the quantile vector or pseudo-random generation vector for the distribution.  Note that the algorithm for the quantile function assumes that the distribution is not bounded from above (i.e., the maximum quantiles is infinity).  In cases where the distribution is bounded from above you would need to vary these algorithms to take account of that additional bound.

Computing the quantile function (vectorised): We are given an input probability vector $\mathbf{p} \in [0,1]^k$ and we want to generate the corresponding vector of quantiles for the distribution.  We assume that we have a function $L$ that can compute the log-CDF of the distribution for any input values.

*

*Define the log-probability vector $\mathbf{l} := \log \mathbf{p}$.


*Define the initial quantile vector $\mathbf{q} := (\infty, ..., \infty)$ (same length as $\mathbf{p}$).


*Define the index set $\mathscr{Q} := \{ i=1,...,k | \mathbf{p}[i] < 1 \}$.


*Define the maximum log-probability for the non-unit probabilities:
$$\ell := \begin{cases}
\max_{i \in \mathscr{Q}} \mathbf{l}[i] & & & \text{if } \mathscr{Q} \text{ is not empty}, \\[6pt]
-\infty & & & \text{if } \mathscr{Q} \text{ is empty}. \\[6pt]
\end{cases}$$

*

*Define the vector of cumulative log-probabilities $\mathbf{L}$ as follows:

*

*Set $r := x_0$, $\mathbf{L} := L(r)$ and $M := L(r)$.

*While $M < \ell$ do:

*

*Set $r := r+1$, $\mathbf{L} := [\mathbf{L}, L(r)]$ and $M := L(r)$.





*For all $i =1,...,k$:

*

*If $\mathbf{l}[i] < 0$ then:

*

*Set $\mathbf{q}[i] := x_0 + \text{sum}(\mathbf{L} < \mathbf{l}[i])$.





*Output the quantile vector $\mathbf{q}$.

Computing the quantile function (vectorised): We are given an input value $k$ and we want to generate a vector of $k$ pseudo-random variables for the distribution.  We assume that we have a function $L$ that can compute the log-CDF of the distribution for any input values.

*

*Generate a vector of pseudo-random values $\mathbf{l} \sim \text{IID log-U}(0,1)$ with length $k$.


*Define the initial quantile vector $\mathbf{q} := (\infty, ..., \infty)$ (with length $k$).


*Define the maximum log-probability $\ell = \max \mathbf{l}$


*Define the vector of cumulative log-probabilities $\mathbf{L}$ as follows:

*

*Set $r := x_0$, $\mathbf{L} := L(r)$ and $M := L(r)$.

*While $M < \ell$ do:

*

*Set $r := r+1$, $\mathbf{L} := [\mathbf{L}, L(r)]$ and $M := L(r)$.





*For all $i =1,...,k$:

*

*If $\mathbf{l}[i] < 0$ then:

*

*Set $\mathbf{q}[i] := x_0 + \text{sum}(\mathbf{L} < \mathbf{l}[i])$.





*Output the quantile vector $\mathbf{q}$.

Programming these algorithms: The above algorithms assume that you have basic programming functionality, including for and while loops and if-then statements.  Depending on the functionality of your programming language, some of the above steps and loops can be "vectorised" in efficient ways.  For example, when programming in R you can make use of various apply functions for certain steps.
