# Approximating a Poisson distribution using a partially observed Gaussian

Note that the problem specification has changed since the original posting. Thanks to whuber for helping me better specify the question: in an attempt to make the question general, I had left out several important constraints.

Problem: Let $y$ be distributed Poisson with rate parameter $\lambda$:

$p(y | \lambda) = \frac{\lambda^y}{y!}e^{-\lambda}$

Let $z^*$ be distributed Gaussian with mean $\mu$ and variance $1$:

$p(z^* | \mu) = \frac{1}{\sqrt{2 \pi}}e^{-\frac{(z^* - \mu)^2}{2}}$

Define $z$ as follows:

$z = k$ if $\gamma_{k} \le z^* < \gamma_{k+1}$, where $\gamma_0 = -\infty$.

Let $\mu = \log\lambda$. How do we choose a single sequence of $\gamma_k$ for $k = 1, 2, \ldots, \infty$, that does not depend on $\lambda$, in order to best approximate the distribution of $y$ with the distribution of $z$, across the full range of possible values of $\lambda \in (0,\infty)$. It is not clear to me exactly what the right criterion is here for operationalizing "best approximating" across the full range of $\lambda$.

Background: The idea here is to do something similar to the Albert and Chib (1993) approach to Gibbs sampling binary and categorical models. Instead of having a single or small number of cutpoints, here we would have an infinite series of cutpoints. The question is how to choose the cutpoints to best approximate a Poisson distribution. As a practical matter, it is pretty easy to draw doubly-truncated normal random variates (Robert 2009).

There is an existing approach to using data augmentation to fit a Poisson regression (Frühwirth-Schnatter and Wagner 2005), but it seems to involve rapidly increasing computation with increasing values of the count variable, which is a very unattractive feature for my application.

I am working on a Gibbs sampler for a hierarchical model to be applied to a very large dataset of count data. I have no particular reason to believe that the distribution of counts are particularly close to Poisson, I am mostly interested in computational efficiency and not having a probability distribution for the counts that is very far from a standard distribution.

References:

James H. Albert and Siddhartha Chib, "Bayesian Analysis of Binary and Polychotomous Response Data"

Sylvia Frühwirth-Schnatter and Helga Wagner, "Data Augmentation and Gibbs Sampling for Regression Models of Small Counts"

Christian Robert, "Simulation of truncated normal variables"

• Please provide full references. The assumption that everyone who reads this recognizes your (name year) references is unlikely to be correct. – Nick Cox Jul 10 '13 at 11:32
• I have added the authors, titles, and links to online copies of the papers. – Ben Lauderdale Jul 10 '13 at 12:05
• Aren't you just asking for the Gaussian cutpoints needed to match the Poisson probabilities? That is, letting $F$ be the CDF of any continuous distribution--which could include any Gaussian you like--you are requiring $F(\gamma_{k+1}) - F(\gamma_k) = p(k|\lambda)$, which immediately gives the full solution inductively. Given the length of your exposition I figure there must be something more going on here, but what is it? – whuber Jul 10 '13 at 14:52
• Maybe that is the answer. I had thought that was the solution, but then I convinced myself that there must be something wrong with that approach. But now I can't remember why I thought that would not do what I wanted. Let me see if I can work out the solution along those lines, I will post it if I can... or clarify the issue if I cannot... – Ben Lauderdale Jul 10 '13 at 15:06
• That sounds hopeless unless you have some constraints on $\lambda$. Even then, you need to stipulate the objective: to minimize the greatest possible deviation? To minimize an expected deviation based on a prior distribution for $\lambda$? The values of $\mu$ and $\sigma$ are irrelevant, by the way, because a solution $(\gamma_k)$ for a standard Gaussian will give the same results as the solution $(\sigma\gamma_k+\mu)$ for the general Gaussian. You might as well just set $\mu=0$ and $\sigma=1$ once and for all. – whuber Jul 10 '13 at 15:40