How to find the covariance of sample mean and sample variance $Cov(M,S^2)$ for Poisson distribution? Suppose that I have a Poisson distribution $P(\lambda)$.
Let $X_1,X_2,\ldots,X_n$ be independent random variables from the distribution mentioned above.
Let us define sample variance $S^2 = \frac{1}{n-1} \sum (X_i - M)^2 $ and sample mean as $M = \frac{1}{n}\sum X_i $.
I want to find the covariance, $Cov(M,S^2)$. I've seen from this answer that when the distribution is symmetric, they are uncorrelated, which makes it zero. I also know that for large $\lambda$ values Poisson distribution is very close to Gaussian distribution, thus becoming symmetric, and probably the covariance is close to zero.
However, I want to find the exact value of $Cov(M,S^2)$, as I am working with small values of $\lambda$

Currently I tried the following:
$Cov(M,S^2) =E( (M-\lambda) (S^2-\lambda) ) = E(MS^2)-\lambda E(M) - \lambda E(S^2) + \lambda^2$
Thus,
$Cov(M,S^2) = E(MS^2) - \lambda^2 - \lambda^2 + \lambda^2  = E(MS^2) -\lambda^2$
I am having trouble with calculating $E(MS^2)$.
Any idea?
 A: One can solve this more generally. For any random variable whose moments exist, the covariance between the sample mean $M$ and the sample variance $S^2$ is:
$$\text{Cov}(M, S^2) = \frac{\mu_3}{n}$$
where $\mu_3$ denotes the $3^{\text{rd}}$ central moment of the parent random variable $X$. 
If $X \sim \text{Poisson}(\lambda)$, then $\mu_3 = \lambda$, so the answer is $\large\frac{\lambda}{n}$. 
A: Let us first center everything, using $\bar X_k=X_k-\lambda$ and $\bar M=M-\lambda$. Then
$$
\mathrm{Cov}(M,S^2)=\mathbb E(\bar MS^2)=\mathbb E(\bar X_1S^2)=\frac1{n-1}\mathbb E(\bar X_1U),
$$
where
$$
U=\sum\limits_{k=1}^n(\bar X_k-\bar M)^2.
$$
Note that $U$ is a linear combination of products $\bar X_k^2$ and $\bar X_k\bar X_i$ for $i\ne k$. Amongst these products, many will not contribute to the expectation of $\bar X_1U$ since $\mathbb E(\bar X_1\bar X_k\bar X_i)=0$ for every $k\ne i$ and $\mathbb E(\bar X_1\bar X_k^2)=0$ for every $k\ne1$. 
Hence, one needs only the coefficient of $\bar X_1^2$ in $U$, which is $c_n=\left(\frac{n-1}n\right)^2+(n-1)\frac1{n^2}=\frac{n-1}n$. This yields $\mathbb E(\bar X_1U)=c_n\mathbb E(\bar X_1^3)$ and $\mathrm{Cov}(M,S^2)=\frac1{n-1}c_n\mathbb E(\bar X_1^3)=\frac1n\mathbb E(\bar X_1^3)$. 
Finally, the third central moment of the Poisson distribution with parameter $\lambda$ is $\mathbb E(\bar X_1^3)=\lambda$ hence
$$
\mathrm{Cov}(M,S^2)=\frac\lambda{n}.
$$
