# Checking Poisson distribution plot using mean and variance relationship?

Im new to Stack exchange, and a student of statistics. Thus I have a question for you guys....

So my questions is asking if it is possible to check the distribution of a random sample comes from a poisson distribution by plotting the mean and the variance of the a distribution.

This is what the code I get for it?

    pos<-rpois(100,5)
x.mean<-rep(0,500)
x.var<-rep(0,500)
for(i in 1:500){
samp<-sample(pos, 10)
x.mean[i]<-mean(samp)
x.var[i]<-var(samp)
}


qqplot-mean and varaince

    quants <- seq(0.01,0.99,length.out=500)
x <- quantile(x.mean, probs=quants)
y <- quantile(x.var, probs=quants)
plot(x,y,xlab="mean",ylab="variance")


And then I notice that the plot(x,y) is pretty much a straight line. So I can conclude that it is a poission distribution. Somehow this just seems wrong though, please let me know what you think?

Best

-
"plot(x,y) is pretty much a straight line. " wouldn't that be also true if you replace the first line by rnorm(100,1,1)? –  user603 Jul 10 '12 at 17:48
When I cut-and-paste the code and run it several times, and compare with @user603's suggestion of a N(1,1) dist'n using the powerful intraocular test (i.e., looking at it), I frankly don't notice much difference between the set of straight lines generated by the Poisson and that generated by the Normal(1,1). The Poisson is perhaps a little straighter, typically, but is also more erratic in the upper tail. –  jbowman Jul 10 '12 at 18:39
@Procrastinator - Only because the OP suggested it; user603's point was, I assume, that the "straight line" is an artifact of the procedure, not a true indicator of Poissionianity, and I was supporting that point. –  jbowman Jul 10 '12 at 18:46
@jb That's right. Moreover, it's not a straight line: throughout much of its length, it is very close to parabolic. –  whuber Jul 10 '12 at 18:49
@Proc You're absolutely right that the sample mean and sample variance will not have the same distribution. Indeed, with a mean of $5$ and samples of $n=10$, the Poisson is already sufficiently close to normal that you can correctly use your experience and intuition about the Normal distribution: the mean will be approximately normal and the variance will be approximately chi-square with $n-1$ df (appropriately scaled). –  whuber Jul 10 '12 at 19:02

There are three issues here:

1. The simulation probably doesn't do exactly what it seems.

2. The analysis does not evaluate whether the relationship between mean and variance is linear.

3. What can one do to check whether data appear Poisson?

Many answers to #3 have appeared elsewhere on our site, so I will focus on the first two issues.

## What the simulation does

The lines of code

pos<-rpois(100,5)
...
samp<-sample(pos, 10)


draw $100$ iid values from a Poisson distribution with mean $5$ and then randomly obtain $10$ of those values (without replacement). Iterating the second command is a form of resampling analysis. It does not really test what happens when data are repeatedly drawn from the same underlying Poisson distribution; instead, it explores the distribution of the $100$ values obtained in the first line, which are fixed for the duration of the simulation.

To study a Poisson distribution, we need to draw repeated independent samples from it:

set.seed(17)      # Create reproducible results
mu <- 3           # Mean
n <- 10           # Sample size
n.trials <- 10000 # Number of trials

sim <- replicate(n.trials, {x <- rpois(n, mu); c(mean(x), var(x))})


At this point sim is a matrix with two rows: one for the means, the next for the variances. Each column summarizes a separate simulation (ergo, there are $10000$ columns in this example).

## What the analysis does

Here it is in a compact form:

xy <- apply(sim, 1, function(x) quantile(x, probs=seq(0.001, 0.999, length.out=99)))


This creates an array of percentiles of the two rows of output, thereby tracing out their distributions. A plot is illuminating:

plot(xy, xlab="mean", ylab="variance",col=hsv(0, 0, 0.25, .4))


Like any Q-Q plot, this one systematically compares percentiles of one distribution (simulation of the variance of samples of $10$ Poisson values) to percentiles of another distribution (simulation of the mean of samples of $10$ Poisson values).

The curve isn't really straight, as we can find by regressing the y values on the x values with a quadratic term included:

z <- xy[,1]^2                  # Quadratic term
fit <- lm(xy[,2] ~ xy[,1] + z) # Regression of y-values against x and x^2
#
# Overplot the fit in red.
x0 <- seq(mu - 2 * sqrt(mu/n), mu + 2 * sqrt(mu/n), length.out=99)
u <- cbind(rep(1,length(x0)), x0, x0^2) %*% coefficients(fit)
lines(x0, u, lwd=2, col=hsv(0, .8, .8, 1))


The quadratic term is highly significant and has a large coefficient:

>coefficients(fit)[3]
z
0.8779665


I encourage interested readers to play with the parameters (mu, n, and n.trials) to see how things (don't) change. Consider also changing rpois in the simulation to something else (such as rnorm).

### Some theory

The mean of $10$ independent draws from many distributions will have approximately a normal distribution. A Poisson (with mean $5$) is no exception: its mean should have a normal distribution of mean $5$ and variance $5/10$.

The CDF of a Poisson distribution with sizable mean (say, $3$ or greater) will itself begin approximating that of a normal distribution, with the approximation improving for larger means. Accordingly, we can hope that the variance of a sample of $10$ independent Poisson draws will begin to assume the distribution of the variance of a sample of $10$ independent draws from a comparable normal distribution. (The variance has a distribution proportional to a chi-squared distribution.) We can check by computing the exact percentiles of these approximating distributions and overplotting them on the simulated values:

q <- (1:100 - 0.5) / 100 # An array of percentiles for plotting
lines(qnorm(q, mean=mu, sd=sqrt(mu/n)), mu*qchisq(q, df=n-1) / (n-1), col="Blue")


Even for this fairly small mean ($5$) and modest sample size ($10$), the fit is extraordinarily good. Unfortunately, this shows that we should expect a similar-looking plot for almost any moderate sample size for many different underlying distributions, thereby demonstrating that this plot tells us little or nothing that is useful concerning whether the data come from a Poisson distribution.

## What can we do instead?

"Checking" whether a distribution has a particular shape can be improved by indicating the ways in which the data might depart from that shape.

For example, Poisson distributions are special examples of negative binomial distributions, which can be considered overdispersed Poisson distributions. Specifically, if we obtain a sequence of Poisson values where the underlying mean of the Poisson distribution changes in a way that approximates a Gamma distribution, then the collected dataset is described by the negative binomial.

In such a setting, we can fit the negative binomial parameters (using maximum likelihood, for instance) and test whether the "non-Poissonness" parameter differs significantly from a value corresponding to a Poisson distribution. In the interest of space, and to keep this reply reasonably elementary, I will not pursue this option here.

Less formally, we can simply apply a Poisson generalized linear model and look at its diagnostics:

set.seed(17)
y <- rnbinom(100, mu=5, size=1)      # NOT Poisson
x <- rpois(100, 5)                   # Poisson
fit.y <- glm(y ~ 1, family=poisson())
fit.x <- glm(x ~ 1, family=poisson())
plot(fit.y)                          # Includes a residual Q-Q plot
plot(fit.x)                          # Includes a residual Q-Q plot


The residual Q-Q plot for the Poisson data is almost linear, whereas that of the negative binomial data has some pronounced extremes

As one can see, it takes a moderate sample size to detect these differences, but at least they do show up.

-
This is great help!! Thank you alot. –  stack_k1 Jul 11 '12 at 15:01