# What is “sampling from a distribution”?

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?

-
It asks of you to generate random numbers so that if you build histogram of large enough sample it would closely match the profile of gamma distribution density, and, moreover, as the size of the sample grows and the bin sizes become smaller the agreement would improve. And yes, you would need a computer for this. – Sasha Oct 6 '11 at 21:47
I suspect it is asking for a description of a method, like this example for selecting a point in a circle, and a demonstration of why the method will work. – Henry Oct 6 '11 at 21:52
I think what is intended is that you give a method, based on "rejection sampling" (this is a standard idea, presumably described there, but you can also search for it), for simulating a gamma distribution. I do not know whether that is supposed to include the grungy details of a computer implementation. – André Nicolas Oct 6 '11 at 21:53
I suspect the point is that you can create a sample from a Gamma distribution with integral $t'$ in various simple ways, such as by adding $t'$ independent (scaled) exponential variates, but there's no obvious way to do this for nonintegral $t$. To accomplish that, generate values from $\Gamma[\lambda', t']$ (where $t'$ is the floor of $t$ and $\lambda'$ is suitably chosen), then reject values by comparing the two PDFs in question. – whuber Oct 6 '11 at 22:35
A solution without proof is here. – Did Oct 7 '11 at 18:42

## 3 Answers

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.

-
I posted recently regarding the Marsaglia Polar Method which is an algorithmic improvement on the Box-Muller Transform. – robjohn Oct 16 '11 at 0:27
@Didier Thanks. – Quinn Culver Oct 16 '11 at 23:52

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.

-

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.

-
I gave this answer a downvote due to being vague (just like the text). Please elaborate. For example, what are the "proper values"? – Quinn Culver Oct 7 '11 at 19:54