# The mode of the Poisson Distribution

Lately, I am doing an investigation on Stirling's formula and its applications. So I thought I could use it to prove that the mode of the Poisson model is approximately equal to the mean. Of course, you do that by considering the curve that is formed by connecting the points of the probabilities of occurrence and the different values of the discrete random variable. Then you differentiate the p.d.f. where for $x!$ you use Stirling's formula $x!\approx \sqrt{2\pi x}~x^xe^{-x}$. The result is $\lnλ-1/(2x)-\ln x$ whose roots cannot be found analytically, but by iterative methods we find that as λ is larger and larger, the mode~mean.

Problem is, I found the following paper online, which seems to be the solution from a Harvard's undergraduate problem set.

It reads "You can also show this by taking the derivative of eq. (2), with Stirling’s expression in place of the $x!$. Furthermore, you can show that $x = a[=λ ~\text{in my case}~]-1/2$ leads to a maximum $P(x)$ value of $P_\max\approx1/\sqrt{2\pi a}$."

Does this puzzle you as much as it puzzles me? My main concern is over the "=" sign: how does this hold? The derivative=0 equation cannot have such an exact solution. Furthermore at such x, how does $P(X=a-1/2)$ give $1/\sqrt{2\pi a}$?

Am I (and my professor) missing something rather obvious or is the solution wrong?

Discuss!

PS: This sort of question might have been asked before, but still, I am really curious that somebody reads the paper in the link above, so that I can figure out what's going on.

-
The link does not work for me. But I note that it contains the word "physics" so it is possible that the solution would strike a mathematician as less than rigorous. (Or as complete nonsense.) – Johan Nov 28 '12 at 14:18
You are right, I had copied and pasted it wrong. It should work now: physics.harvard.edu/academics/undergrad/probweek/sol84.pdf – Ryuky Nov 28 '12 at 14:20

The weight $w(n)$ of the Poisson distribution with positive parameter $\lambda$ at the integer $n\geqslant0$ is $w(n)=\mathrm e^{-\lambda}\lambda^n/n!$ hence $w(n+1)/w(n)=\lambda/(n+1)$. One sees that $w(n)\gt w(n-1)$ for every $n\lt\lambda$ while $w(n+1)\lt w(n)$ for every $n\gt\lambda-1$. Thus, the mode of the Poisson distribution with parameter $\lambda$ is the highest integer $n_\lambda$ such that $n\lt\lambda$.

The estimation of $w_\lambda=\max\limits_nw(n)$ when $\lambda\to\infty$ is direct through Stirling's equivalent since $\lambda-1\lt n_\lambda\lt \lambda$, and indeed yields $\lim\limits_{\lambda\to\infty}\sqrt{2\pi\lambda}\cdot w_\lambda=1$.

Edit: Extend the sequence $(w(n))_{n\geqslant0}$ to a function $W$ defined on $\mathbb R^+$ through the formula $W(x)=\mathrm e^{-\lambda}\lambda^x/\Gamma(x+1)$. Thus, $W(n)=w(n)$ for every nonnegative integer $n$, and the function $W$ is maximal at $x_\lambda$ such that $\log\lambda=\psi(x_\lambda+1)$, where $\psi$ denotes the polygamma function. The asymptotic expansion of $\psi$ at infinity is $\psi(z)=\log(z)-\frac1{2z}+o(\frac1z)$ hence $\mathrm e^{\psi(z)}=z-\frac12+o(1)$ and $\lambda=x_\lambda+\frac12+o(1)$, that is, $\lim\limits_{\lambda\to\infty}x_\lambda-\lambda=-\frac12$.

-
To rely on Stirling's approximation to compute the mode would be, as somebody put it on this page, complete nonsense. – Did Nov 28 '12 at 14:27
I understand your reasoning about the mode. However, did you read the paper from the link? The part I pointed out is false, right? – Ryuky Nov 28 '12 at 14:32
Can you please elaborate on the last bit? – Ryuky Nov 28 '12 at 14:47
See Edit.   – Did Nov 28 '12 at 14:59
@Did are you able to comment on the possibility of a solution for the mode in the bivariate case? stats.stackexchange.com/q/133083/7258 – Killercam Jan 16 '15 at 14:24

To find the mode of the Poisson distribution, for $k > 0$, consider the ratio $$\frac{P\{X = k\}}{P\{X = k-1\}} = \frac{e^{-\lambda}\frac{\lambda^k}{k!}}{e^{-\lambda}\frac{\lambda^{k-1}}{(k-1)!}} = \frac{\lambda}{k}$$ which is larger than $1$ for $k < \lambda$ and smaller than $1$ for $k > \lambda$.

• If $\lambda < 1$, then $P\{X = 0\} > P\{X = 1\} > P\{X > 2\} \cdots$ and so the mode is $0$.

• If $\lambda > 1$ is not an integer, then the mode is $\lfloor\lambda\rfloor$ since $P\{X = \lceil\lambda\rceil\} < P\{X = \lfloor\lambda\rfloor\}$.

• If $\lambda$ is an integer $m$, then $P\{X = m\} = P\{X = m-1\}$ and so either $m$ or $m-1$ can be taken to be the mode.

In all cases, the mode and the mean differ by less than $1$. You do not need to use Stirling's approximation at all.

-

If all you're trying to prove is that the mode of the Poisson distribution is approximately equal to the mean, then bringing in Stirling's formula is swatting a fly with a pile driver. You have $$P(X=x) = \frac{\lambda^x e^{-\lambda}}{x!}.$$ The mean is $\lambda$. Now let us seek the mode.

Observe that $$\frac{P(X=x+1)}{P(X=x)} = \frac{\lambda^{x+1}/(x+1)!}{\lambda^x/x!} = \frac{\lambda}{x+1}.$$ This is $\ge 1$ if $x\le\lambda-1$ and $\le1$ if $x\ge \lambda-1$. Thus $$P(X=x+1)\quad \left.\begin{cases} \ge \\ = \\ \le \end{cases}\right\}\quad P(X=x)$$ according as $$x\quad \left.\begin{cases} \le \\ = \\ \ge \end{cases}\right\}\quad \lambda-1.$$ The mode is therefore the integer part of $\lambda-1$. Except when $\lambda$ is an integer, in which case two consecutive integers are both modes.

The linked item at Harvard was trying to do much more than just find the mode.

Later addendum: Although the mode of the distribution must be within the set that is the support of the distribution, which is $\{0,1,2,3,\ldots\}$, the linked paper seeks the value of $x$ that maximizes $\lambda^2 e^{-\lambda}/x!$ when non-integer values of $x$ are allowed. We could define $x!=\Gamma(x+1)$. However, the author uses Stirling's approximation, and that at least allows us to do something in closed form. Consider (at this point I'll call it $a$ instead of $\lambda$) $$\frac{a^x e^{-a}}{x!} \approx \frac{a^x e^{-a}}{x^x e^{-x}\sqrt{2\pi x}}.$$ Since $e^{-a}$ and $\sqrt{2\pi}$ don't depend on $x$ we can lose those and look at $$a^x e^x x^{-x} x^{-1/2} = \exp\left((x\log a) + x - x\log x - \frac12 \log x\right).$$ Since $\exp$ is an increasing function, we can seek the value of $x$ that maximizes the expression inside it, and that will be the value of $x$ that makes the derivative of that expression $0$: $$\frac{d}{dx} \left( (x\log a) + x - x\log x - \frac12 \log x \right) = \log a - \log x -\frac{1}{2x} = \log\left(\frac a x\right) - \frac{1}{2a}\left(\frac a x\right) = 0.$$ I don't see a way to get this in closed form, unless you want to bring in Lambert's $W$. That being the case, I don't know why he didn't just use $x!=\Gamma(x+1)$. But I entered this command into Wolfram Alpha:

f(x) = log(6) - log(x) - 1/(2x); x from 4 to 6

Lo and behold: it crosses the $x$-axis very close to $5.5$.

So the paper is not very explicit, to say the least, about how this conclusion was arrived at.

Still later addendum: Now I've entered this command into Wolfram Alpha:

f(x) = 6^x*e^(-6) / Gamma(x+1); from 5.49 to 5.51

It looks as if the maximum is near $5.494$.

-
Okay, I agree with all that, but if I understand correctly, in the part that I have quoted, he is saying that the maximum value of $P$ occurs at $x=a-1/2$. How is this so? – Ryuky Nov 28 '12 at 14:41
Yes. But except for that, how did he come up with this value of $x$. He is saying that he is taking the derivative, which after setting equal to $0$ gives a rather complicated equation that cannot be solved exactly. So, please tell me, what is he talking about? And if you plug this value of $x$ into the Poisson how can you get $1/\sqrt{2\pi a}$ There might be something that I am missing since I am still at high school, or what he is doing does not make sense. – Ryuky Nov 28 '12 at 15:00
@Ryuky : I've added some material on how $a-1/2$ was arrived at. The answer could be: numerically. But that's only for particular values of $a$. – Michael Hardy Nov 28 '12 at 17:26
The answer could be: numerically... No, this is a consequence of the expansion of the polygamma function $\psi$, see my answer. – Did Nov 28 '12 at 18:12
A minor correction: if $\lambda$ is not an integer, the mode is the integer part of $\lambda$ (not of $\lambda-1$ as you have it), cf. Didier's answer (and mine). – Dilip Sarwate Nov 28 '12 at 20:28