5
$\begingroup$

Exercise 4.11.3 of Grimmett and Stirzaker's Probability and Random Processes reads "Use the rejection method to sample from the gamma density $\Gamma(\lambda,t)$ where $t (\geq 1)$ may not be assumed integral."

What exactly is this exercise asking for? More generally, what is meant by "sampling from a distribution" (or "sampling from a density")? Section 4.11 of the given textbook refers to "sampling" over and over, but doesn't define it.

It seems, based on the examples, that it means to simply find a random variable with the given distribution. If that's correct, why is it called "sampling"?

Am I supposed to use a computer for this exercise?

  • 0
    Have you looked in the **[master class](http://cg.scs.carleton.ca/~luc/rn$b$ooki$n$dex.html)**? (**Hint**: IX.) Regarding bounties, last I $k$$n$ew, the author mai$n$tains an open bounty of some Belgian beer for an algorithm to generate gamma random variates in [one line of code](http://cg.scs.carleton.ca/~luc/wsc96.ps). (Unfortunately, for whatever reason, I can't seem to find a reference to that bounty at the moment.)2011-10-15

3 Answers 3

4

In view of the context, no computer is involved. Rather, one asks for a method to generate a random variable $X$ with a given distribution (in your case, Gamma) from a collection, possibly infinite, of independent random variables, which are either uniformly distributed over $(0,1)$, or some simple transformations of these.

Assume for example that one asks for a standard gaussian random variable. A popular approach, called the Box-Muller method, is to use $U_1$ and $U_2$ i.i.d. uniform on $(0,1)$ and to set $X=\sqrt{-2\log U_1}\cos(2\pi U_2)$. See here for more explanations and examples.

In the Box-Muller method, one uses $U_1$ and $U_2$ to get one random variable $X$. Acception-rejection methods are different because they use a random number of random variables to get $X$. Their general principle is to use a given number of i.i.d. random variables, say two uniform random variables $U_1$ and $U_2$, and to test if these satisfy a well chosen property, say $\mathcal Q(U_1,U_2)$. If $\mathcal Q(U_1,U_2)$ holds, one accepts $(U_1,U_2)$ in the sense that one sets $X=\Phi(U_1,U_2)$ for a suitable function $\Phi$ and the game is over. Otherwise, that is, if $\mathcal Q(U_1,U_2)$ does not hold, one rejects $U_1$ and $U_2$, and one performs the same test with two other uniform random variables $U_3$ and $U_4$, say. If $\mathcal Q(U_3,U_4)$ holds, one accepts $(U_3,U_4)$ in the sense that one sets $X=\Phi(U_3,U_4)$ and the game is over. Otherwise one rejects $U_3$ and $U_4$ and one turns to $U_5$ and $U_6$, and so on. One can also keep some information from the rejected random values, for example their number, to make the value $X$, but the procedure above is always the rough idea.

Hence, in rejection methods, the total number of random variables $U_n$ used to produce one value of $X$ is random. In your case, you are given a possibly infinite sequence $(U_n)_{n\geqslant1}$ of i.i.d. uniform random variables and you are asked to devise an algorithm to produce one random variable $X$ with the desired Gamma distribution, using almost surely a finite number of random variables $U_n$.

A good and relatively short introduction is the first chapter (The acceptance rejection method with applications; generating a standard normal random variable) of the lecture notes of this course. Another source, written for high energy astrophysicists and remarkably down-to-earth, is here.

The second link provides the code of an algorithm to generate Gamma distributions from what they call a Lorentzian distribution, which is also known as a standard Cauchy distribution. This is quite convenient since the latter is simply the distribution of $\tan(\pi U)$ with $U$ uniform on $(-1,1)$. I would rather leave you the pleasure of putting these bits and pieces together to concoct your algorithm, unless you get stuck at some point.

  • 0
    I [posted recently](http://math.stackexchange.com/questions/69245/transform-uniform-distribution-to-normal-distribution-using-lindeberglevy-clt/69284#69284) regarding the Marsaglia Polar Method which is an algorithmic improvement on the Box-Muller Transform.2011-10-16
2

The question is asking you to devise a scheme to generate random numbers according to the specified distribution, Gamma$(\lambda, t)$. Formally, we are coming up with a way to construct a random variable (i.e. a measurable function from a probability space to the reals) with a predetermined distribution given only a sequence of independent random variables of predetermined distributions, using a particular method - in this case, so called "rejection sampling". The usual rejection sampling algorithm to get samples from a distribution with density $f$, using another density $g$, is along the lines of this:

  1. Sample $Y \sim g(y)$ from some candidate distribution, and $U \sim \mbox{Uniform}(0, 1)$ independently.
  2. Accept $Y = X$ as a simulation from $f(x)$ if $U \le \frac{f(Y)}{Mg(Y)}$ where $M \ge \sup_x \frac{f(x)}{g(x)}$.
  3. Otherwise, go back to step 1, using random variables which are independent of those which have already been generated.

More precisely, we have sequences $Y_1, Y_2, ...$ and $U_1, U_2, ...$ of independent random variables on a probability space $(\Omega, \mathcal F, P)$, and if you undertake this algorithm for each fixed $\omega \in \Omega$ then the random variable $X$ can be shown to be distributed according to $f$, provided that the support of $g$ contains the support of $f$ and that such an $M < \infty$ exists. As Dider notes, we want our algorithm to terminate after a finite number of steps almost surely, and it can be shown that this one does.

In the case of the Gamma distribution, you can take $g(y)$ to be another Gamma density that you know how to sample from, say Gamma$(\lambda', t')$ where t' is integral. You can easily get draws from such a Gamma using only uniform random variables, and then turn them into draws from the Gamma$(\lambda, t)$ using the above method.

-3

I believe that you would need to use something like MATLAB for that. You need to input the proper values and you can generate a set of random values that you can work from.

  • 0
    I gave this answer a downvote due to being vague (just like the text). Please elaborate. For example, what are the "proper values"?2011-10-07