Peak of Poisson distribution I was experimenting with Poisson distribution on desmos and noticed that as $\lambda$ (the mean of distribution) gets larger, the x coordinate of distribution peak approaches $\lambda-0.5$. Why is it so?
 A: The mode of the Poisson distribution occurs at the value $\text{Mode}(\lambda) = \lceil \lambda \rceil - 1$, whereas your proposed approximation is:
$$\widehat{\text{Mode}}(\lambda) \equiv \lambda - \tfrac{1}{2}.$$
Letting $u(\lambda) \equiv \lceil \lambda \rceil - \lambda$ denote the "upper remainder" of the number $\lambda$ we can write the difference between the true mode and your approximation as:
$$\begin{align}
\text{DIFF} \equiv
\text{Mode}(\lambda) - \widehat{\text{Mode}}(\lambda)
&= (\lceil \lambda \rceil - 1) - (\lambda-\tfrac{1}{2}) \\[6pt]
&= (\lceil \lambda \rceil - \lambda) -\tfrac{1}{2} \\[6pt]
&= u(\lambda) - \tfrac{1}{2}. \\[6pt]
\end{align}$$
It is simple to establish that $-\tfrac{1}{2} \leqslant \text{DIFF} < \tfrac{1}{2}$ so your approximation is near to the true mode.  As $\lambda \rightarrow \infty$ the relative error in the approximation converges to zero.
A: The analysis at https://stats.stackexchange.com/a/211612/919 shows the mode of any Poisson$(\lambda)$ distribution is near $\lambda$ itself.  Although that question concerns only integral $\lambda,$ its results answer the present question, too.
Let $p_\lambda(k)$ be the Poisson$(\lambda)$ probability for $k\in\{0,1,2,\ldots,\}$ (all other probabilities are zero, of course).  These probabilities are proportional to the ratios
$$p_\lambda(k) \ \propto\  \frac{\lambda^k}{k!}.$$
As pointed out in the foregoing link, this means two successive Poisson probabilities are related by
$$p_\lambda(k+1) = \frac{\lambda}{k+1} p_\lambda(k).$$
Consequently, as $k$ progresses from $0$ through all integers less than $\lambda-1,$ the probability increases from $p_\lambda(k)$ to $p_\lambda(k+1);$ and then once $k$ exceeds $\lambda-1,$ the probability decreases.  Therefore

All Poisson probability functions $p_\lambda$ rise to a peak and then fall again.  The peak (a mode) occurs at the greatest integer less than $\lambda,$ written $\lfloor \lambda\rfloor.$  When $\lambda$ is an integer, the peak occurs at the two neighboring values $\lfloor \lambda\rfloor$ and $\lfloor \lambda\rfloor + 1.$

To illustrate, here is a plot of the modes based on a brute-force search.

The search was conducted by this R function.  It evaluates the Poisson probabilities for $k$ between two suitable extreme quantiles, finds the indexes where the largest probability occurs, and returns an array of the values of $k$ corresponding to those indexes.
mode <- function(lambda) {
  q <- min(1/2, ppois(floor(lambda), lambda))          # This is *a* probability
  k <- do.call(seq, as.list(qpois(c(q, 1-q), lambda))) # Search limits
  p <- dpois(k, lambda)
  k[which(p == max(p))]
}

