ktp.tex
\documentclass[11pt]{article}
\usepackage{amsmath}
\usepackage{indentfirst}
\DeclareMathOperator{\var}{var}
\newcommand{\Prmu}{\Pr\nolimits_\mu}
\begin{document}
\title{The $K$-Truncated Poisson Distribution}
\author{Charles J. Geyer}
\maketitle
\section{Introduction}
This document works through the details of the $k$-truncated Poisson distribution,
a special case of which is the zero-truncated Poisson distribution. The $k$-truncated Poisson
distribution is the distribution of a Poisson random variable $Y$ conditional on the event
$Y > k$. It has one parameter, which we may take to be $\mu = E(Y)$.
Since $\mu$ is not the mean (or anything else simple) of the distribution of $Y$ conditioned
on the even $Y > k$, we do not call $\mu$ the mean, rather we call it
the \emph{original parameter}.
If $f_\mu$ is the probability mass function (PMF) of $Y$, then the
PMF $g_\mu$ of the $k$-truncated Poisson distribution is
defined by
\begin{equation} \label{eq:pmf-one}
g_\mu(x)
=
\frac{f_\mu(x)}{1 - \sum_{j = 0}^k f_\mu(j)},
\qquad x = k + 1, k + 2, \ldots
\end{equation}
Plugging in the formula for the Poisson PMF, we get
\begin{equation} \label{eq:pmf-two}
\begin{split}
g_\mu(x)
& =
\frac{\frac{\mu^x}{x !} e^{- \mu}}{1 - \sum_{j = 0}^k \frac{\mu^j}{j !} e^{- \mu}}
\\
& =
\frac{\mu^x}{x ! (e^\mu - \sum_{j = 0}^k \frac{\mu^j}{j !})}
\end{split}
\end{equation}
\section{Exponential Family Properties}
Clearly \eqref{eq:pmf-one} is the PMF of a one parameter exponential family having
canonical statistic $x$ and canonical parameter $\theta = \log(\mu)$.
Of course, the original parameter is $\mu = \exp(\theta)$.
The cumulant function for the family is then
\begin{equation} \label{eq:psi}
\psi(\theta)
=
\log \left(e^{e^\theta} - \sum_{j = 0}^k \frac{e^{j \theta}}{j !} \right)
\end{equation}
and has derivatives
\begin{equation} \label{eq:tau}
\begin{split}
\tau(\theta)
=
\psi'(\theta)
& =
\frac{ e^\theta \cdot e^{e^\theta}
- \sum_{j = 1}^k \frac{e^{j \theta}}{(j - 1) !} }
{e^{e^\theta} - \sum_{j = 0}^k \frac{e^{j \theta}}{j !}}
\\
& =
\frac{ \mu
- e^{- \mu} \sum_{j = 1}^k \frac{\mu^j}{(j - 1) !} }
{1 - e^{- \mu} \sum_{j = 0}^k \frac{\mu^j}{j !}}
\end{split}
\end{equation}
and
\begin{equation} \label{eq:psi-double-prime}
\begin{split}
\psi''(\theta)
& =
\frac{ (e^\theta + e^{2 \theta}) \cdot e^{e^\theta}
- \sum_{j = 1}^k \frac{j e^{j \theta}}{(j - 1) !} }
{e^{e^\theta} - \sum_{j = 0}^k \frac{e^{j \theta}}{j !}}
-
\tau(\theta)^2
\\
& =
\frac{ (\mu + \mu^2) - e^{- \mu}
\sum_{j = 1}^k \frac{j \mu^j}{(j - 1) !} }
{1 - e^{- \mu} - \sum_{j = 0}^k \frac{\mu^j}{j !}}
-
\tau(\theta)^2
\end{split}
\end{equation}
By exponential family theory we know $\psi'(\theta) = E_\theta(X)$ and
$\psi''(\theta) = \var_\theta(X)$, where $X$ is the canonical statistic.
Thus from our definition of $\tau(\theta)$
in \eqref{eq:tau} it follows that $\tau'(\theta) = \psi''(\theta) > 0$ for all $\theta$.
Hence the map $\tau$ is one-to-one and defines an invertible change of parameter.
Since $\tau(\theta) = E_\theta(X)$, it is called the \emph{mean value parameter}.
It is the mean of the distribution under discussion: $k$-truncated Poisson.
\subsection{Check}
If these are correct, then \eqref{eq:tau} should be $E(X)$ and
the fraction in \eqref{eq:psi-double-prime} should be $E(X^2)$ when $X$ has
the $k$-truncated Poisson distribution.
\begin{align*}
E(X)
& =
\sum_{x = k + 1}^\infty x g_\mu(x)
\\
& =
\sum_{x = k + 1}^\infty \frac{x f_\mu(x)}{1 - \sum_{j = 0}^k f_\mu(j)}
\\
& =
\frac{\mu - \sum_{j = 0}^k j f_\mu(j)}{1 - \sum_{j = 0}^k f_\mu(j)}
\end{align*}
agrees with \eqref{eq:tau}.
\begin{align*}
E(X^2)
& =
\sum_{x = k + 1}^\infty \frac{x^2 f_\mu(x)}{1 - \sum_{j = 0}^k f_\mu(j)}
\\
& =
\frac{\mu + \mu^2 - \sum_{j = 0}^k j^2 f_\mu(j)}{1 - \sum_{j = 0}^k f_\mu(j)}
\end{align*}
agrees with \eqref{eq:psi-double-prime}.
\subsection{Computing}
As always, we wish to compute things, in this case the cumulant function and its
first two derivatives, without overflow or cancellation error.
Problems arise when $\mu$ is nearly zero or when $\mu$ is very large.
\subsubsection{Cumulant Function}
From \eqref{eq:psi} we get, using $\mu = \exp(\theta)$,
\begin{equation} \label{eq:psi-comp}
\begin{split}
\psi(\theta)
& =
\mu + \log \left(1 - e^{- \mu} \sum_{j = 0}^k \frac{\mu^j}{j !} \right)
\\
& =
\mu + \log \Prmu \{ Y > k \}
\end{split}
\end{equation}
where $Y \sim \text{Poi}(\mu)$. This looks fairly stable whether $\mu$ is large or small.
We will leave the calculation of the log Poisson probability to R.
\subsubsection{First Derivative of Cumulant Function}
From \eqref{eq:tau} we get, using $\mu = \exp(\theta)$,
\begin{subequations}
\begin{equation} \label{eq:tau-comp}
\begin{split}
\tau(\theta)
& =
\frac{ \mu - e^{- \mu} \sum_{j = 1}^k \frac{\mu^j}{(j - 1) !} }
{\Prmu \{ Y > k \}}
\\
& =
\frac{ \mu \left[
1 - e^{- \mu} \sum_{j = 0}^{k - 1} \frac{\mu^j}{j !} \right] }
{\Prmu \{ Y > k \}}
\\
& =
\frac{ \mu \Prmu \{ Y \ge k \}}
{\Prmu \{ Y > k \}}
\\
& =
\mu \left( 1 + \frac{\Prmu \{ Y = k \}}
{\Prmu \{ Y > k \}} \right)
\end{split}
\end{equation}
While this looks good for large $\mu$ it is not at all clear that it behaves well when
$\mu$ is small. As $\mu \to 0$ (and $\theta \to - \infty$) we
have $\tau(\theta) \to k + 1$. Let us see if we can get a computationally stable way
to compute that without using L'Hospital's rule.
\begin{equation} \label{eq:tau-comp-too}
\begin{split}
\tau(\theta)
& =
\mu
+
\frac{\mu \Prmu \{ Y = k \}}
{\Prmu \{ Y > k \}}
\\
& =
\mu
+
\frac{\mu^{k + 1} e^{- \mu} / k !}
{\mu^{k + 1} e^{- \mu} / (k + 1) ! + \Prmu \{ Y > k + 1 \}}
\\
& =
\mu
+
\frac{\displaystyle k + 1}
{\displaystyle 1 +
\frac{\Prmu \{ Y > k + 1 \}}{\Prmu \{ Y = k + 1 \}}}
\end{split}
\end{equation}
\end{subequations}
When $\mu$ is nearly zero, then the fraction in the denominator is also nearly zero
and we get nearly $k + 1$ with no chance of overflow. Oops! It can produce \verb@NaN@
(IEEE not a number) when the fraction in the denominator is $0 / 0$,
actually underflow over underflow. If we special case this, then everything works.
Actually, our second formula \eqref{eq:tau-comp-too}, seems too work just as well
as \eqref{eq:tau-comp}, even when $\mu$ is very large. In calculating $\tau(\theta)$
from zero to 1000 in steps of 0.1 both formulas give the same answers to within machine
precision (relative error about $10^{- 16}$) whenever they do not give \verb@Inf@,
which they do for precisely the same arguments $\theta \le 709.7$.
In hindsight, this is no surprise. When $\mu \simeq \infty$, then
$\Prmu \{ Y > k + 1 \} \approx 1$ and $\Prmu \{ Y = k + 1 \} \approx 0$
and the quotient in the denominator of \eqref{eq:tau-comp-too} either is very large or overflows
giving \verb@Inf@ when IEEE arithmetic is in use (what happens on ancient computers without it
is problematic), and the whole fraction is nearly zero. Hence, when $\mu \simeq \infty$,
\eqref{eq:tau-comp-too} adds something very small or zero to $\mu$.
\subsubsection{Second Derivative of Cumulant Function}
We start our computation of $\psi''(\theta)$ by noting that $\psi''(\theta) = \tau'(\theta)$,
and, because $d \mu / d \theta = \mu$,
$$
\tau'(\theta) = \mu \frac{\tau(\mu)}{d \mu}.
$$
Thus we differentiate our ``good'' expression \eqref{eq:tau-comp-too} for $\tau$
expressed in terms of $\mu$. It will simplify notation if we also define
$$
\beta
=
\frac{\Prmu \{ Y > k + 1 \}}{\Prmu \{ Y = k + 1 \}}
=
\frac{e^\mu \Prmu \{ Y > k + 1 \}}{e^\mu \Prmu \{ Y = k + 1 \}}
$$
and note that \eqref{eq:tau-comp-too} says
$$
\tau = \mu + \frac{k + 1}{1 + \beta}.
$$
so
$$
\frac{d \tau}{d \mu} = 1 - \frac{k + 1}{(1 + \beta)^2} \cdot \frac{d \beta}{d \mu}.
$$
To calculate $d \beta / d \mu$ we first figure out
$$
\frac{d}{d \mu} e^\mu \Prmu \{ Y > k + 1 \}
=
\frac{d}{d \mu} \sum_{y = k + 2}^\infty \frac{\mu^y}{y !}
=
\sum_{y = k + 2}^\infty \frac{\mu^{y - 1}}{(y - 1) !}
=
e^\mu \Prmu \{ Y > k \}
$$
and
$$
\frac{d}{d \mu} e^\mu \Prmu \{ Y = k + 1 \}
=
\frac{d}{d \mu} \frac{\mu^{k + 1}}{(k + 1) !}
=
\frac{\mu^k}{k !}
=
e^\mu \Prmu \{ Y = k \}
$$
Then
\begin{align*}
\frac{d \beta}{d \mu}
& =
\frac{d}{d \mu} \frac{e^\mu \Prmu \{ Y > k + 1 \}}{e^\mu \Prmu \{ Y = k + 1 \}}
\\
& =
\frac{e^\mu \Prmu \{ Y > k \}}{e^\mu \Prmu \{ Y = k + 1 \}}
-
\frac{e^\mu \Prmu \{ Y > k + 1 \}}{\left( e^\mu \Prmu \{ Y = k + 1 \} \right)^2}
\cdot
e^\mu \Prmu \{ Y = k \}
\\
& =
\frac{\Prmu \{ Y > k \}}{\Prmu \{ Y = k + 1 \}}
-
\frac{\Prmu \{ Y > k + 1 \}}{\Prmu \{ Y = k + 1 \}}
\cdot
\frac{\Prmu \{ Y = k \}}{\Prmu \{ Y = k + 1 \}}
\\
& =
\frac{\Prmu \{ Y > k + 1 \}}{\Prmu \{ Y = k + 1 \}} + 1
-
\frac{\Prmu \{ Y > k + 1 \}}{\Prmu \{ Y = k + 1 \}}
\cdot
\frac{\Prmu \{ Y = k \}}{\Prmu \{ Y = k + 1 \}}
\\
& =
\beta + 1
-
\beta \cdot
\frac{\Prmu \{ Y = k \} }{\Prmu \{ Y = k + 1 \}}
\\
& =
\beta + 1 - \beta \cdot \frac{k + 1}{\mu}
\end{align*}
So finally
\begin{equation} \label{eq:psi-double-prime-comp}
\psi''(\theta)
=
\mu \left( 1 - \frac{k + 1}{1 + \beta}
\left( 1 - \frac{k + 1}{\mu} \cdot \frac{\beta}{1 + \beta} \right) \right)
\end{equation}
\section{Random Variate Generation}
To simulate a $k$-truncated Poisson distribution, the simplest method is to simulate
ordinary Poisson random variates (using the \verb@rpois@ function in R) and reject all
of the simulations less than or equal to $k$.
This works well unless $\mu = \exp(\theta)$, the mean of the untruncated Poisson distribution
is nearly zero, in which case the acceptance rate is also nearly zero.
In that case, another simple rejection sampling scheme, simulates $Y \sim \text{Poi}(\mu)$
and uses $X = Y + m$ as the rejection sampling proposal, where $m$ is a nonnegative integer
(the case $m = 0$ is the case already discussed).
The ratio of target density to proposal density is
\begin{equation} \label{eq:rat}
\begin{split}
\frac{g_\mu(x)}{f_\mu(y)}
& =
\frac{g_\mu(x)}{f_\mu(x - m)}
\\
& =
\frac{f_\mu(x) I(x > k)}{ f_\mu(x - m) \left( 1 - \sum_{j = 0}^k f_\mu(j) \right)}
\\
& =
\frac{(x - m) ! \cdot \mu^m I(x > k)}{ x ! \left( 1 - \sum_{j = 0}^k f_\mu(j) \right)}
\end{split}
\end{equation}
where $I(x > k)$ is one when $x > k$ and zero otherwise.
This achieves its upper bound (considered as a function of $x$) when $x = \max(m, k + 1)$.
To avoid the ``max'' let us impose the condition that $m \le k + 1$, so the max is achieved
when $x = k + 1$ and is
\begin{equation} \label{eq:maxrat}
\frac{(k + 1 - m) ! \cdot \mu^m}{ (k + 1) ! \left( 1 - \sum_{j = 0}^k f_\mu(j) \right)}
\end{equation}
Thus we accept proposals with probability \eqref{eq:rat} divided by \eqref{eq:maxrat},
which is
\begin{equation} \label{eq:accept-prob}
\frac{(x - m) ! (k + 1) !}{ x ! (k + 1 - m) !} \cdot I(x > k)
\end{equation}
As noted at the beginning of the discussion, when $m = 0$ we accept a proposal $x$ with
probability $I(x > k)$.
When $m = 1$ we accept a proposal $x$ with probability
$$
\frac{k + 1}{x} \cdot I(x > k)
$$
and so forth.
To understand the performance of this algorithm, hence to understand how to chose $m$, we need
to calculate the acceptance rate
\begin{align*}
a(k, m)
& =
E \left\{
\frac{(X - m) ! (k + 1) !}{ X ! (k + 1 - m) !} \cdot I(X > k)
\right\}
\\
& =
\frac{(k + 1) !}{ (k + 1 - m) !} \cdot
E \left\{
\frac{Y !}{ (Y + m) !} \cdot I(Y > k - m)
\right\}
\\
& =
\frac{(k + 1) !}{ (k + 1 - m) !}
\sum_{y = k + 1 - m}^\infty
\frac{y !}{ (y + m) !} \cdot \frac{\mu^y}{y !} e^{-\mu}
\\
& =
\frac{(k + 1) !}{ (k + 1 - m) !}
\cdot \frac{1}{\mu^m}
\sum_{y = k + 1 - m}^\infty
\frac{\mu^{y + m}}{(y + m) !} e^{-\mu}
\\
& =
\frac{(k + 1) !}{ (k + 1 - m) !}
\cdot \frac{1}{\mu^m}
\sum_{w = k + 1}^\infty
\frac{\mu^w}{w !} e^{-\mu}
\\
& =
\frac{(k + 1) !}{ (k + 1 - m) !}
\cdot \frac{1}{\mu^m}
\cdot \Pr(Y > k)
\end{align*}
Everything is fixed in our formula for acceptance rate except $m$ which we many choose
to be any integer $0 \le m \le k + 1$. Consider
$$
\frac{a(k, m + 1)}{a(k, m)}
=
\frac{(k + 1 - m)}{\mu}
$$
this is greater than one (so it pays to increase $m$) when
$$
k + 1 - m < \mu
$$
which suggests we make
$$
m = \lceil k + 1 - \mu \rceil
$$
so long as this denotes a nonnegative integer (otherwise we set $m = 0$).
The performance of this algorithm seems to be fine for small $k$.
However the worst case acceptance rate, which occurs for $\mu$ between $k / 4$ and $k / 2$,
does seem to go to zero as $k$ goes to infinity.
For a zero-truncated Poisson distribution the worst case acceptance rate is 63.2\%.
For a two-truncated Poisson distribution the worst case acceptance rate is 48.2\%.
For a twenty-truncated Poisson distribution the worst case acceptance rate is 21.6\%.
For a one-hundred-truncated Poisson distribution the worst case acceptance rate is 10.2\%.
\end{document}
