Asymptotic Moments of the Binomial Distribution, $E(X/(np))^k = 1 + O(k^2/n)$?

The $k$th moment of a Binomial$(n,p)$ random variable $X$ is $$\mathbb{E}(X^k)=\sum_{j=1}^k{k\brace j}j!{n\choose j}p^j,\tag1$$ where ${k\brace j}$ is the Stirling number of the second kind.

For large $n$, we have $j!{n\choose j}\approx n^j-{j\choose 2}n^{j-1}$. Taking the two leading terms in (1) gives the approximation $$\mathbb{E}(X^k)\approx \left[n^k-{k\choose 2}n^{k-1}\right]p^k+{k\choose 2}n^{k-1}p^{k-1}=(np)^k+{k\choose 2}(np)^{k-1}(1-p),$$ since ${k\brace k}=1$ and ${k\brace k-1}={k\choose 2}.$

The right universal bound is $$E\left[\left(\frac{X}{np}\right)^k\right] \le \left(1+\frac{k}{2np}\right)^k.$$ For $k = O(\sqrt{np})$ this explains the $1+O(k^2/(np))$ behaviour.

Using the Moment Generating Function

Let $$m(t) = E[e^{t X}] = (1-p+pe^t)^n \le \exp(\mu(e^t-1)),$$ where $\mu=np$. This last bound is the moment generation of a Poisson random variable, so the following bound will hold in that case as well.

Now $E[X^k] \le m(t)t^{-k}k!$ follows easily from the expansion $E[e^{tX}]=\sum_i E[X^i]t^i/i!$. However, we will need the slightly stronger bound: $$E[X^k] \le m(t)(k/(et))^k.$$ This follows from the basic inequality $1+x\le e^x$, where we substitute $tx/k-1$ for $x$ to get $tx/k \le e^{tx/k-1} \implies x^k \le e^{tx}(k/(et))^k$. Taking expectations we get the intended bound.

We define $B=k/\mu$ like above, and take $t$ such that $t e^t=B$. (This $t$ is also known as $W(B)$, using the Lambert function.)

We then have $$\begin{align} E[(X/\mu)^k] &\le m(t)(\mu t)^{-k}(k/e)^k \\&\le \exp(\mu(e^t-1))\big(\frac{k}{e\mu t}\big)^k \\&= \exp(\mu(B/t-1)+tk-k) \\&= \exp(k f(B)) , \end{align}$$ where $f(B) = 1/t+t-1-1/B$. We can bound $\exp(f(x))$ by $B/\log(B+1) \le 1+B/2 \le e^{B/2}$ using standard bounds on the Lambert function.

This finally gives the bounds $$ \|X/\mu\|_k \le \frac{B}{\log(B+1)} \le e^{B/2}. $$ Taking $k$ powers this means $$ \text{and}\quad E[(X/\mu)^k] \le \exp(k^2/(2np)), $$ where $e^{k^2/(2np)} = 1+O(k^2/(np))$ when $k = O(\sqrt{np})$.

Old attempt using norms:

I found the following interesting non-asymptotic bound:

Let $\|X\|_k = E[|X|^k]^{1/k}$ and use the triangle inequality $$\|X\|_k \le np + \|X-np\|_k.$$

To bound the second term, we use Latala's formula for sums of iid. random variables: $$\|X\|_k = \|\sum_i X_i - E[X_i]\|_k \lesssim \inf_{\max\{2,\frac kn\}\le s\le k} \frac ks \left(\frac nk\right)^{1/s} \|X_i - E[X_i]\|_s$$

For Bernoulli random variables with probability $p$: $\|B-p\|_k = \left(p (1-p)^k+(1-p) p^k\right)^{1/k} \le p^{1/k}$.

Let $\beta = \frac{k}{np}$. We assume $k\le n$ and $np\ge 1$. Optimizing in Latala's we have

$$ \|X\|_k \lesssim \begin{cases} \sqrt{n p k} & \text{if $\beta < e$} \\ \frac{k}{\log\beta} & \text{if $\beta\ge e$}. \end{cases} $$

Putting it all together we have

$$ E\left[\left(\frac{X}{np}\right)^k\right] \le \begin{cases} (1 + C\sqrt{\frac{k}{np}})^k & \text{if $\beta < e$} \\ (1 + C\frac{k}{n p\log\beta})^k & \text{if $\beta\ge e$}. \end{cases} $$

In particular for $k = O((n p)^{1/3})$ we get $E(\frac{X}{np})^k=1 + O(k^3/(n p))$. This is not quite as good as the conjectured bound, but at least it is rigorous and works even when $p$ is also allowed to depend on $n$.

It also shows that the moments of $\frac{X}{np}-1$ are within a constant of the equivalent Gaussian, as long as $k < np$.

Likely one has to avoid the triangle inequality to show the stronger conjecture.