3
$\begingroup$

In Stochastic Simulation: Algorithms and Analysis by Søren Asmussen, on Page 38

A Poisson r.v. $N$ with rate $\lambda$ ($P(X = n) = e^{−\lambda} \frac{\lambda^n}{n!}$) can be constructed using the relation between the exponential distribution and the Poisson process: generate $X_1 , X_2 , ...$ as i.i.d. exponentials with rate $\lambda$ one at a time until the random time at which the sum exceeds 1, say $X_1 + · · ·+ X_N < 1 < X_1 + · · ·+ X_{N +1}$ . Then N is the desired Poisson r.v. To reduce the number of evaluations of special functions, use $N = max \lbrace n \geq 0 : \prod_{i=1}^n U_i > e^{-\lambda} \rbrace$

where $U_i$ is r.v. uniform in $[0,1]$.

I was wondering why to "reduce the number of evaluations of special functions", the $N$ in the last equation has Poisson distribution with rate λ?

Thanks and regards!

  • 0
    "reduce the number of evaluations of special functions" - because it's not *cheap* to evaluate special functions. Even the logarithm will take a while to evaluate if you have to do it a number of times, which is par for the course for a rejection method.2010-09-24

1 Answers 1

5

Note that by inverse sampling we can generate exponential r.v. with rate $\lambda$ by

$X=F^{-1}(U)$ where $F(x) = 1-e^{-\lambda x}$, is the c.d.f of exponential r.v. with rate $\lambda$ and $U\sim Uni[0,1]$. Thus $F^{-1}(U)=\frac{-\log(1-U)}{\lambda}$.

But $1-U\sim Uni[0,1]$ so if $U_i\sim Uni[0,1]$ then

$X_i=\frac{-\log(U_i)}{\lambda}$ is exponential with rate $\lambda$.

$X_1+\ldots +X_n= -\frac{1}{\lambda}(\log(U_1)+\ldots +\log(U_n)) = -\frac{1}{\lambda}\log(U_1U_2\ldots U_n)$

That is $X_1+\ldots + X_n<1$ iff $\prod_{i=1}^nU_i>e^{\lambda}$. Clearly in both cases $N$ is actually the same random variable.

Why reduce the evaulation of special functions?

Suppose $\lambda=10^{-3}$ (small). In the first method, where you generate each $X_i$ in turn add them up and check if totals less than one, you would expect to generate $10^3$ exponentials. Each one requires you to evaluate the logarithm (with argument in $(0,1)$) and simulate a uniform r.v.. So you have 1 thousand evaluations of a special function, in the region where it is hard to approximate.

In the second method, you expect to generate $10^3$ uniform r.v. and check if the product is greater than $e^{-10^3}$. This method requires you to calculate $e^{-10^3}$ but only once.

Both methods generate what you desired but the second method will be much more efficient on a computer.