In short, the function I am interested in is the following:
$$I_n(p) = \sum_{k=0}^n \binom{n}{k} p^k (1-p)^{n-k} \left[h(p) - h\left(\tfrac{k}{n}\right)\right],$$
where $h(x) \triangleq -x \log x - (1-x) \log (1-x)$ is the standard binary entropy function. I am able to prove that for large $n$ this function attains a local maximum around $p \approx \frac{1.338}{n}$ at which $I_n(p) \approx \frac{0.84}{n}$ in case binary logarithms are used. By assuming that $p = \frac{\alpha}{n}$ for $\alpha = O(1)$ and using a Poisson approximation of the binomial distribution, the value $\alpha \approx 1.338$ is a (numerically found) maximizing value of
$$\max_{\alpha > 0} \left\{\sum_{z=1}^{\infty} \frac{\alpha^z}{z!} e^{-\alpha} (z \log z) - \alpha \log \alpha\right\}.$$
Note that the latter expression does not contain $n$ anymore; I am interested in large-$n$ asymptotics of $I_n(p)$ which leads to the above, $n$-independent expression. 

What I would like to show is that this point $p = \frac{\alpha}{n}$ is actually not just a local, but a global maximum, i.e., there are no values $p = \omega(\frac{1}{n})$ with even higher asymptotic values of $I_n(p)$ for large $n$. So my question can be formulated as:

> Can we prove that for large $n$, $\max\limits_{p \in [0,1]} I_n(p) \to I_n(\frac{\alpha}{n})$ with $\alpha \approx 1.338184$? 

Numerical inspection of $I_n(p)$ for say $n = 100\,000$ shows that $I_n(p) \propto \frac{1}{n}$ is nice and smooth, has two global maxima at $\frac{\alpha}{n}$ and $1 - \frac{\alpha}{n}$ with value $\frac{0.84}{n}$ and a local minimum at $p = \frac{1}{2}$ with value $\frac{0.72}{n}$ (and two trivial minima at $0$ and $1$ with value $0$). The figure below sketches what $I_n(p)$ (scaled with a factor $n$) looks like for small $n$. As $n$ increases, the function might be best described as a flat line at height $0.72$ with two peaks close to $0$ and $1$ of height $0.84$.

![enter image description here][1]

I've tried various approaches using e.g. Taylor expansions of the entropy, but I am not able to prove that "order terms" are actually small for arbitrary $p$. Also, a bound given [here][2] using Jensen's inequality is not sharp enough; it's roughly a factor $2$ too loose to prove the above statement.

Any help would be appreciated!

---

## Attempt 1 ##
For large $n$, the dominating terms in the summation in $I_n$ come from values $k \approx np$. Expanding $h(\frac{k}{n})$ around $k = np$ leads to
$$h\left(\frac{k}{n}\right) = h(p) + \left(\frac{k}{n} - p\right) h'(p) + \frac{1}{2}\left(\frac{k}{n} - p\right)^2 h''(p) + \frac{1}{6}\left(\frac{k}{n} - p\right)^3 h^{(3)}(p) + \dots.$$
Substituting this back into the definition of $I_n(p)$, this leads to
$$I = \sum_{k=0}^n \binom{n}{k} p^k (1-p)^{n-k} \left[-\left(\tfrac{k}{n} - p\right) h'(p) - \tfrac{1}{2}\left(\tfrac{k}{n} - p\right)^2 h''(p) - \tfrac{1}{6}\left(\tfrac{k}{n} - p\right)^3 h^{(3)}(p) - \dots\right].$$
Now the term $\frac{k}{n}$ leads to a factor $\frac{np}{n} = p$ when pulled out of the summation, so the first term disappears. For the second term, using the first and second moments of the binomial distribution, the term $\frac{k^2}{n^2}$ turns into a $p(1-p+np)/n = p(1-p)/n + p^2$, the cross-term $-2 k p/n$ turns into $-2p^2$, and the term $p^2$ stays $p^2$. So in total, the second term results in $\frac{1}{2n}  p(1-p)h''(p)$. As $h'(p) = \log_2(\frac{1-p}{p})$ and $h''(p) = \frac{-1}{p(1-p)\ln 2}$, this term in the end contributes $\frac{1}{2n \ln 2} \approx \frac{0.72}{n}$ as expected. So we have:
$$I_n(p) = \frac{1}{2n \ln 2} + \sum_{k=0}^n \binom{n}{k} p^k (1-p)^{n-k} \left[- \tfrac{1}{6}\left(\tfrac{k}{n} - p\right)^3 h^{(3)}(p) - \dots\right].$$
At this point, what remains is to show that all remaining order terms are small, if $np(1-p) = \omega(1)$. 

  [1]: https://i.sstatic.net/Imyoj.png
  [2]: http://web.media.mit.edu/~bcroy/papers/bcroy_expected-entropy.pdf