Normal approximation to the Poisson distribution Here in Wikipedia it says:

For sufficiently large values of $λ$, (say $λ>1000$), the normal distribution with mean $λ$ and variance $λ$ (standard deviation $\sqrt{\lambda}$), is an excellent approximation to the Poisson distribution. If $λ$ is greater than about 10, then the normal distribution is a good approximation if an appropriate continuity correction is performed, i.e., $P(X ≤ x),$ where (lower-case) $x$ is a non-negative integer, is replaced by $P(X ≤ x + 0.5).$
$F_\mathrm{Poisson}(x;\lambda) \approx F_\mathrm{normal}(x;\mu=\lambda,\sigma^2=\lambda)$

Unfortunately this isn't cited. I want to be able to show / prove this with some rigour. How can you actually say the normal distribution is a good approximation when $\lambda > 1000$, how do you quantify this 'excellent' approximation, what measures were used?
The furthest I've got with this is here where John talks about using the Berry–Esseen theorem and approximates the error in the two CDFs. From what I can see he does not try any values of $\lambda \geq 1000$.
 A: Suppose $X$ is Poisson with parameter $\lambda$, and $Y$ is normal with mean and variance $\lambda$.  It seems to me that the appropriate comparison is between $\Pr(X = n)$ and $\Pr(Y \in [n-\frac12,n+\frac12])$.  Here for simplicity I write $n = \lambda  + \alpha \sqrt\lambda$, that is, we are interested when $n$ corresponds to $\alpha$ standard deviations from the mean.
So I cheated.  I used Mathematica.  So both $\Pr(X = n)$ and $\Pr(Y \in [n-\frac12,n+\frac12])$ are asymptotic to
$$ \frac 1{\sqrt{2\pi \lambda}} e^{-\alpha^2/2} $$
as $\lambda \to \infty$.  But their difference is asymptotic to
$$ \frac{\alpha  \left(\alpha ^2-3\right) e^{-\alpha ^2/{2}}}{6 \sqrt{2
   \pi } \lambda } $$
If you plot this as a function of $\alpha$, you will get the same curve as is shown in the second to last figure in http://www.johndcook.com/blog/normal_approx_to_poisson/.
Here are the commands I used:
  n = lambda + alpha Sqrt[lambda];
  p1 = Exp[-lambda] lambda^n/n!;
  p2 = Integrate[1/Sqrt[2 Pi]/Sqrt[lambda] Exp[-(x-lambda)^2/2/lambda], {x, n-1/2, n+1/2}];
  Series[p1, {lambda, Infinity, 1}]
  Series[p2, {lambda, Infinity, 1}]

Also, with a bit of experimentation, it seems to me that a better asymptotic approximation to $\Pr(X = n)$ is $\Pr(Y \in [n-\alpha^2/6,n+1-\alpha^2/6])$.  Then the error is
$$ -\frac{\left(5 \alpha ^4-9 \alpha ^2-6\right) e^{-{\alpha ^2}/{2}}
   }{72 \sqrt{2 \pi } \lambda ^{3/2} } $$
which is about $\sqrt\lambda$ times smaller.
A: The derivation from the binomial distribution might gain you some insight. 
We have a binomial random variable; 
$$ p(x) = {n \choose x} p^x (1-p)^{n-x}$$
This can alternatively be computed recursively; 
$$ p(x) = \frac{(n-x+1)p}{x(1-p)}p(x-1)$$ 
If you keep the initial condition;
$$ p(0) = (1-p)^n $$
Now let us assume that $n$ is large and $p$ is small but the average success of $p(x)$ is constant $(np = \lambda)$. Then we can do the following;
$$ P( X = i ) = {n \choose i} p^x (1-p)^{n-x} $$
We use that $p = \lambda / n$. 
$$ P( X = i ) = \frac{n!}{(n-i)!i!} \left(\frac{\lambda}{n}\right)^i \left(1-\frac{\lambda}{n}\right)^{n-i}  $$
We switch some variables around and evaluate;
$$  P( X = i ) = \frac{n(n-1)(n-2)\cdots(n-i+1)}{n^i} \frac{\lambda^i}{i!} \frac{(1-\frac{\lambda}{n})^n}{(1-\frac{\lambda}{n})^i} $$ 
From calculus we know that $ \lim_{n\to\infty} (1 + x/n)^n = e^x $. We also know that $[n(n-1)(n-2)\cdots(n-i+1)]/n^i \approx 1 $ because both the top and bottom are polynomials of degree $i$. 
This leads to the conclusion that as $ n \to \infty$: 
$$ P(X=i) \to \frac{ e^{-\lambda}{\lambda^i}}{i!} $$ 
You can then verify that $E(X) = \lambda$ and $\operatorname{Var}(X) = \lambda$ via the definition. We know that the binomial distribution approximates the normal under the conditions of the De Moivre-Laplace Theorem as long as you correct for the continuity, which is why $P(X\le x)$ is replaced by $P(X\le x+0.5)$.
A: Glen_b is correct in that "good fit" is a very subjective notion. However, if you want to verify that your poisson distribution is reasonably normal, you can use a hypothetical Kolmorgov-Smirnov test with the null hypothesis being $H_{0}:$ The CDF came from a $N(\lambda,\lambda)$ distribution, assuming your sample will come from a poisson($\lambda$). Since you are not actually testing a sample, but one distribution against another, you need to think carefully about the sample size and significance level you assume for this hypothetial test (since we are not using the KS test in its typical fashion). That is:


*

*Pick a representative, hypothetical sample size, n, and adjust the significance level of the test to a typical value, e.g., 5% . 


Now, calculate the Type II error rate for this test assuming your data actually come from a poisson($\lambda$). Your degree of fit with a normal distribution will be this Type II error rate, in the sense that samples of size n from your particular poisson distribution will, on average, be accepted $\beta$% of the time by a KS normality test at your selected significance level. 
Anyway, that's just one way to go about getting a sense of "goodness of fit". However, all rely on some subjective notions of "goodness" that you will have to define for yourself.
