Past exam of the simulation and statistical programming course of the University of Oxford / 2021 / 2

For the exponential distribution with rate $\lambda$, the cumulative distribution function is $F(x)=1-e^{-\lambda x}$ for $x\ge 0$, and zero for $x<0$. Its quantile function isFor inverse transform sampling, draw $U$ from the uniform distribution on $(0,1)$ and return $-\log (1-U)/\lambda$. Equivalently return $-\log U/\lambda$, since $1-U$ has the same uniform distribution. Indeed, for $x\ge 0$, $\mathbb{P}(-\log U/\lambda >x)=\mathbb{P}(U<e^{-\lambda x})=e^{-\lambda x}$.

Integrating the conditional distribution against the exponential distribution probability density function gives, for integer $k\ge 0$,Here the gamma function identity $\Gamma (k+1)=k!$ evaluates the integral. Henceon $\{0,1,2,\dots \}$. Put $r=1/(1+\lambda )$ and $\rho =1-r$. The geometric series gives $F(k)=1-r^{k+1}$. Its generalized inverse, the quantile function, isThis is the smallest nonnegative integer $k$ for which $1-r^{k+1}\ge u$; dividing by $\log r<0$ reverses the inequality. Draw $U$ from the uniform distribution on $(0,1)$ and use this formula. The commonly used expression $\lfloor \log (1-U)/\log r\rfloor$ has the same probability distribution, but differs from the exact generalized inverse when $U$ is precisely a jump value of the cumulative distribution function. Such values have probability zero for an ideal continuous uniform draw.

Write $t=1-\rho \in (0,1)$. The target Poisson distribution and proposal geometric distribution have probability mass functions $p_k=e^{-\lambda }{\lambda }^k/k!$ and $q_k=(1-t)t^k$. Their ratio isThe successive ratio $R_{k+1}/R_k=\lambda /[t(k+1)]$ shows that a maximizing integer is $m=\lfloor \lambda /t\rfloor$, with a tie at $m-1$ when $\lambda /t$ is an integer. Thus the tight envelope isFor any fixed $t$, draw $Y\sim \mathrm{Geom}(1-t)$ using (ii)'s inversion formula with this new parameter, and draw an independent $U$ from the uniform distribution on $(0,1)$. Accept whenotherwise restart. The acceptance probability is $1/M(t)$ and the expected number of proposals is $M(t)$.

To minimize that expected count, set $K=\lceil \lambda \rceil$. The exact answer isHere is a proof that includes the changes of maximizing integer. Each $\log R_k(t)$ is a convex function of $t$, since its second derivative is $(1-t{)}^{-2}+k/t^2>0$. Their pointwise maximum $\log M(t)$ is also a convex function. When $m$ is fixed its derivative isIf $K-1<\lambda \le K^2/(K+1)$, take $t_{\ast }=\lambda /K$, where $m$ changes from $K$ on the left to $K-1$ on the right. The left derivative is nonpositive because $\lambda \le K^2/(K+1)$; the right derivative is nonnegative because $\lambda \ge K-1$. This is a minimum. If $K^2/(K+1)<\lambda \le K$, take $t_{\ast }=K/(K+1)$. There $K<\lambda /t_{\ast }\le K+1$ and the active branch has $m=K$ except for the endpoint tie; its derivative vanishes. Convexity again proves the global minimum. The formulas agree at the boundary.

For a fully specified optimal algorithm, compute $K,t_{\ast },{\rho }_{\ast }$ and $m_{\ast }=\lfloor \lambda /t_{\ast }\rfloor$ once, then use the preceding acceptance test with $t=t_{\ast }$ and $m=m_{\ast }$. At an integer tie either maximizing integer gives the same test. In code, factorials can be evaluated through the logarithm of the gamma function:and acceptance is $\log U\le \min (0,\log a_Y)$. This avoids overflowing factorials.

The optimum minimizes expected trials, the usual idealization of time when each proposal and test has constant cost. Actual runtime also depends on the implementation of discrete inversion and factorial evaluation. Simply minimizing an upper bound for $M(t)$ would not necessarily give the exact optimum above; see optimal geometric envelope for a Poisson distribution.

Because $U\in [0,1]$, the conditional acceptance probability is $a_k=\min (1,{\lambda }^k)$, even when the printed threshold exceeds one. Repeating rejected trials gives the clipped acceptance reweighting identity:For $0<\lambda <1$, $a_k={\lambda }^k$, including $a_0=1$. The denominator is $e^{-\lambda }{\sum }_k{\lambda }^{2k}/k!=e^{{\lambda }^2-\lambda }$, soFor $\lambda \ge 1$ every proposal is accepted, so reweighting by ${\lambda }^k$ would be incorrect. At $\lambda =1$ both formulas agree. Termination has probability one since the acceptance probability is $e^{{\lambda }^2-\lambda }>0$ in the first case and one in the second.

Define $a_y=\min \{1,{\lambda }^y/7\}$. Given $Y=y$, the direct-return branch has probability $a_y$, while the fallback branch distributes its remaining mass uniformly over $0,\dots ,y$. The law of total probability therefore gives, for integer $x\ge 0$,The summation includes $y=x$: a failed direct-return test can still return $x$ through the discrete uniform distribution draw. There is no rejection loop and no additional normalization. Summing the displayed probability mass function over $x$ gives ${\sum }_y{p}_y[a_y+(1-a_y)(y+1)/(y+1)]=1$. At $\lambda =1$, for example, $a_y=1/7$ at every $y$; it does not become one. This is an instance of the accept-or-fallback mixture identity.