# Using the rejection sampling with the method of inversion

I am hoping to write some rejection algorithm code in R to approximate a $\text{Gamma}(k,\lambda)$ distribution.

The problem is more for educational purposes than real-world implementation.

Given an $\text{Exponential}(\lambda)$ distribution with PDF:

$$f_{\lambda}(x) = \lambda\text{exp}(-\lambda x)$$

the CDF can be expressed as:

$$F_{\lambda}(x) = 1 - \text{exp}(-\lambda x)$$

where $x>0$ and $\lambda>0$.

And the inverse of the CDF can be expressed as:

$$F^{-1}(U) = -\text{log}(1-U)/\lambda$$

where $U$ is a random number generated such that: $U \sim \text{U}[0,1]$.

Assuming that:

$$X_{1}, X_{2}, \dots, X_{k} \stackrel{\text{ iid }}{ \sim }\text{Exponential}(\lambda)$$

and:

$$Y = X_{1} + X_{2} + \dots + X_{k}$$

then:

$$Y \sim \text{Gamma}(k, \lambda)$$

Now, the above approach for simulating a $\text{Gamma}(k,\lambda)$ from a sum of $k$ $\text{Exponential}(\lambda)$ random variables will, naturally, only work where $k \in {\bf N}$, but let's assume the objective is to simulate $\text{Gamma}(k,\lambda)$ where $k \notin {\bf N}$. It has been suggested to use an accept-reject approach with an envelope function being the $\text{Gamma}(\lfloor{k}\rfloor,\lambda-1)$ density, where $\lfloor{k}\rfloor$ is the function floor(k).

EDIT:

In the above example, it's clear that target distribution is given by:

$$f(x) = \frac{\lambda^{k}}{\Gamma(k)}x^{k-1}e^{-\lambda x}$$

But in order to solidify these concepts, could somebody provide the form of the simulating distribution ($h(x)$), the envelope function ($g(x) = M*h(x)$) and, of course, the optimal value of $M$?

Solidifying these ideas would really help me out.

### Closely related

Rejection sampling for a modified Gamma distribution

How to quickly sample X if exp(X) ~ Gamma? (one answer provides a clever R implementation).

-
I am afraid there is some confusion there about the nature of the accept-reject method. For one thing, you cannot use an exponential to approximate Gamma distributions with $\alpha<1$. For another, the idea of using a sum of $n$ exponential variates will produce a Gamma$(n,\lambda)$ variate. I suggest you look at our Example 2.19 (page 52) in our Monte Carlo Statistical Methods book with George Casella. –  Xi'an Oct 24 '12 at 20:50
(followed from above) Precisely, I do not understand your remark about $f(x)<x\sim G(\alpha,\beta)$. This is not the standard condition: $u<f(x)/Mg(x)$. Inversion is an altogether different method: do you want to use it for simulating the exponential? It is trivial: $X=-\log(U)/\lambda\sim E(\lambda)$ is de facto an inversion. –  Xi'an Oct 25 '12 at 9:15
@Xi'an, thanks for your input. My apologies. As you can probably tell, this area is brand new to me, and, originally, the question was very poorly constructed; since then, I have revised the question. I hope it makes more sense now. –  user9171 Oct 25 '12 at 13:08
This is the example we detail in our book, indeed. –  Xi'an Oct 25 '12 at 13:22
@Xi'an, fantastic! I'll be sure to check it out! It must be a common example. Thanks again! –  user9171 Oct 25 '12 at 13:24
Given the target density $$f(x) = \frac{\lambda^{k}}{\Gamma(k)}x^{k-1}e^{-\lambda x} \quad k>1$$ let us write $k_0=\lfloor{k}\rfloor$ and take as envelope density $$g(x) = \frac{\lambda_0^{k_0}}{\Gamma(k_0)}x^{k_0-1}e^{-\lambda_0 x}$$ with $\lambda_0>0$. Then the ratio $f/g$ is given by $$\frac{\lambda^{k}}{\Gamma(k)}\,\frac{\Gamma(k_0)}{\lambda_0^{k_0}}\, x^{k-k_0}\,e^{-(\lambda -\lambda_0)x}$$ which is bounded in 0 because $k-k_0\ge 0$ and at $+\infty$ when $\lambda -\lambda_0>0$. Under this assumption, the maximum is obtained at $x^*$ solution of $$\frac{k-k_0}{x} = \lambda -\lambda_0$$ by taking the derivative of the ratio. Therefore $$x^* = \frac{k-k_0}{\lambda -\lambda_0}$$ and \begin{align*} M &= f(x^*)/g(x^*)\\ &= \frac{\lambda^{k}}{\Gamma(k)}\,\frac{\Gamma(k_0)}{\lambda_0^{k_0}}\, (x^*)^{k-k_0}\,e^{-(\lambda -\lambda_0)x^*}\\ &= \frac{\lambda^{k}}{\Gamma(k)}\,\frac{\Gamma(k_0)}{\lambda_0^{k_0}}\, (x^*)^{k-k_0}\,e^{-(k-k_0)} \end{align*} i.e. $$M = \frac{\lambda^{k}}{\Gamma(k)}\,\frac{\Gamma(k_0)}{\lambda_0^{k_0}}\, (x^*/e)^{k-k_0}$$ As described in Example 2.19 (page 52) in our Monte Carlo Statistical Methods book with George Casella, you can further optimise the choice of $\lambda_0$ by minimising $M$ in $\lambda_0$, which leads to $\lambda_0/k_0=\lambda=k$, i.e. $$\lambda_0=\frac{k_0\lambda}{k}$$