Okay, last time I said that the next post would be about logarithmic Sobolev inequalities… well, it’s not. I lied. Instead, we’re taking a momentary break from concentration inequalities. It’s time for some relaxation of the convex variety.

Actually, that’s a bit of lie too because the post isn’t about convex relaxation, but it is about convex optimization. Andiamo!

Maximum likelihood of Poisson distributions

We’re going to work through a problem from Convex Optimization by Boyd and Vandenberghe. (If you’re curious it’s problem 7.7.) Here’s the set up.

• We have $X_1,\dots,X_n$ independent Poisson random variables. That is each $X_i$ has a probability mass function parameterized by $\lambda_i > 0$ that is given by

$$\prob (X_i = k) = \frac{e^{-\lambda_i}\lambda_i^k}{k!}.$$

• The data $x_1,\dots,x_n$ represent the number of times that one of $n$ possible events occurred (independently of one another).
• We have $m$ detectors and event $i$ is recorded by detector $j$ with probability $p_{ji}$. We assume that $p_{ji} \ge 0$, $\sum_{j=1}^m p_{ji} \le 1$, and that the $p_{ji}$ are given.
• The total number of events recorded by detector $j$ is

$$y_j = \sum_{i=1}^n y_{ji}, \quad j = 1,\dots,m.$$

• Our goal is to estimate the $\lambda_i$ based on observations $y_j$ and the given probabilities $p_{ji}$ by solving a convex optimization problem.

To pose this estimation problem as a convex optimization problem, we’ll use the fact that the $y_{ji}$ were generated from a Poisson distribution with mean $p_{ji}\lambda_i$. Intuitively, this means we only capture a fraction, in particular $p_{ji}$, of the Poisson-distributed events associated with $\lambda_i$. For those who have had a course on Markov chains this should cause your thinned Poisson process neurons to fire.

Now, we use the fact that the sum of Poisson random variables is itself a Poisson random variable. In particular, we have that the $y_j$ were generated by a Poisson distribution with parameter

$$p_j^T \lambda = \sum_{i=1}^n p_{ji}\lambda_i,$$

where $p_j \in \reals^n$ and $\lambda \in \reals_+^n$ is the vector of $\lambda_i$.

The likelihood of the data has the form

$$\prod_{j=1}^m \frac{(p_j^T \lambda)^{y_j}e^{-(p_j^T \lambda)}}{y_j!}$$

which gives a log-likelihood of

$$\sum_{j=1}^m \left( y_j\log(p_j^T \lambda)-p_j^T \lambda-\log(y_j!) \right).$$

We’ll throw away the $-\log(y_j!)$ term because it is a constant. So, the estimation problem is the convex optimization problem

$$\begin{array}[ll] \mbox{\text{maximize}} & \sum_{j=1}^m \big( y_j \log(p_j^T \lambda) - p_j^T \lambda \big) \end{array}$$

with optimization variable $\lambda \in \reals_+^n$ and problem data $y_1,\dots,y_m$ and $p_{ji}$ for $i = 1,\dots,n$, $j = 1,d\dots,m$. More explicitly, we have

$$\begin{array}[ll] \mbox{\text{maximize}} & \sum_{j=1}^m \big( y_j \log(p_j^T \lambda) - p_j^T \lambda \big) \\ \mbox{subject to} & \lambda \succeq 0. \end{array}$$

Put it to use. Now, go out and find yourself a particle accelerator, equip it with some detectors, start colliding particles and count your Higgs-Bosons and make some inferences.

Next up

… logarithmic Sobolev inequalities—for real this time (well, next time)!

References

This post is based on material from…

• Convex Optimization by Stephen Boyd and Lieven Vandenberghe. Check it out here: Convex Optimization.