4
$\begingroup$

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).

$\endgroup$
10
  • 2
    $\begingroup$ 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. $\endgroup$
    – Xi'an
    Oct 24, 2012 at 20:50
  • $\begingroup$ (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. $\endgroup$
    – Xi'an
    Oct 25, 2012 at 9:15
  • $\begingroup$ @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. $\endgroup$
    – user9171
    Oct 25, 2012 at 13:08
  • $\begingroup$ This is the example we detail in our book, indeed. $\endgroup$
    – Xi'an
    Oct 25, 2012 at 13:22
  • $\begingroup$ @Xi'an, fantastic! I'll be sure to check it out! It must be a common example. Thanks again! $\endgroup$
    – user9171
    Oct 25, 2012 at 13:24

1 Answer 1

4
$\begingroup$

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} $$

$\endgroup$
1
  • $\begingroup$ Do you mean $\lambda_0/k_0=\lambda/k$ near the end there, rather than $\lambda_0/k_0=\lambda=k$? $\endgroup$
    – Glen_b
    Nov 19, 2019 at 23:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.