What is the distribution of the rounded down average of Poisson random variables? If I have random variables $X_1,X_2,\ldots,X_n$ that are Poisson distributed with parameters $\lambda_1, \lambda_2,\ldots, \lambda_n$, what is the distribution of $Y=\left\lfloor\frac{\sum_{i=1}^n X_i}{n}\right\rfloor$ (i.e. the integer floor of the average)?
A sum of Poissons is also Poisson, but I am not confident enough in statistics to determine if it is the same for the case above.
 A: Y will not be Poisson. Note that Poisson random variables take on non negative integer values.  Once you divide by a constant you create a random variable that can have non-integer values.  It will still have the shape of the Poisson.  It is just that the discrete probabilities may occur at non-integer points.
A: A generalization of the question asks for the distribution of $Y = \lfloor X/m \rfloor$ when the distribution of $X$ is known and supported on the natural numbers.  (In the question, $X$ has a Poisson distribution of parameter $\lambda = \lambda_1 + \lambda_2 + \cdots + \lambda_n$ and $m=n$.)
The distribution of $Y$ is easily determined by the distribution of $mY$, whose probability generating function (pgf) can be determined in terms of the pgf of $X$.  Here's an outline of the derivation.

Write $p(x) = p_0 + p_1 x + \cdots + p_n x^n + \cdots$ for the pgf of $X$, where (by definition) $p_n = \Pr(X=n)$.  $mY$ is constructed from $X$ in such a way that its pgf, $q$, is
$$\eqalign{q(x) &=& \left(p_0 + p_1 + \cdots + p_{m-1}\right) + \left(p_m + p_{m+1} + \cdots + p_{2m-1}\right)x^m + \cdots + \\&&\left(p_{nm} + p_{nm+1} + \cdots + p_{(n+1)m-1}\right)x^{nm} + \cdots.}$$
Because this converges absolutely for $|x| \le 1$, we can rearrange the terms into a sum of pieces of the form
$$D_{m,t}p(x) = p_t + p_{t+m}x^m + \cdots + p_{t + nm}x^{nm} + \cdots$$
for $t=0, 1, \ldots, m-1$.  The power series of the functions $x^t D_{m,t}p$ consist of every $m^\text{th}$ term of the series of $p$ starting with the $t^\text{th}$: this is sometimes called a decimation of $p$.  Google searches presently don't turn up much useful information on decimations, so for completeness, here's a derivation of a formula.
Let $\omega$ be any primitive $m^\text{th}$ root of unity; for instance, take $\omega = \exp(2 i \pi / m)$.  Then it follows from $\omega^m=1$ and $\sum_{j=0}^{m-1}\omega^j = 0$ that
$$x^t D_{m,t}p(x) = \frac{1}{m}\sum_{j=0}^{m-1} \omega^{t j} p(x/\omega^j).$$
To see this, note that the operator $x^t D_{m,t}$ is linear, so it suffices to check the formula on the basis $\{1, x, x^2, \ldots, x^n, \ldots \}$.  Applying the right hand side to $x^n$ gives
$$x^t D_{m,t}[x^n] =  \frac{1}{m}\sum_{j=0}^{m-1} \omega^{t j} x^n \omega^{-nj}=  \frac{x^n}{m}\sum_{j=0}^{m-1} \omega^{(t-n) j.}$$
When $t$ and $n$ differ by a multiple of $m$, each term in the sum equals $1$ and we obtain $x^n$. Otherwise, the terms cycle through powers of $\omega^{t-n}$ and these sum to zero.  Whence this operator preserves all powers of $x$ congruent to $t$ modulo $m$ and kills all the others: it is precisely the desired projection.
A formula for $q$ follows readily by changing the order of summation and recognizing one of the sums as geometric, thereby writing it in closed form:
$$\eqalign{
q(x) &= \sum_{t=0}^{m-1} (D_{m,t}[p])(x) \\
     &= \sum_{t=0}^{m-1} x^{-t} \frac{1}{m} \sum_{j=0}^{m-1} \omega^{t j} p(\omega^{-j}x ) \\
     &= \frac{1}{m} \sum_{j=0}^{m-1} p(\omega^{-j}x) \sum_{t=0}^{m-1} \left(\omega^j/x\right)^t  \\
     &= \frac{x(1-x^{-m})}{m} \sum_{j=0}^{m-1} \frac{p(\omega^{-j}x)}{x-\omega^j}.
}$$
For example, the pgf of a Poisson distribution of parameter $\lambda$ is $p(x) = \exp(\lambda(x-1))$.  With $m=2$, $\omega=-1$ and the pgf of $2Y$ will be
$$\eqalign{
q(x) &= \frac{x(1-x^{-2})}{2} \sum_{j=0}^{2-1} \frac{p((-1)^{-j}x)}{x-(-1)^j} \\
     &= \frac{x-1/x}{2} \left(\frac{\exp(\lambda(x-1))}{x-1} + \frac{\exp(\lambda(-x-1))}{x+1}\right) \\
     &= \exp(-\lambda) \left(\frac{\sinh (\lambda  x)}{x}+\cosh (\lambda  x)\right).
}$$
One use of this approach is to compute moments of $X$ and $mY$.  The value of the $k^\text{th}$ derivative of the pgf evaluated at $x=1$ is the $k^\text{th}$ factorial moment.  The $k^\text{th}$ moment is a linear combination of the first $k$ factorial moments.  Using these observations we find, for instance, that for a Poisson distributed $X$, its mean (which is the first factorial moment) equals $\lambda$, the mean of $2\lfloor(X/2)\rfloor$ equals $\lambda- \frac{1}{2} + \frac{1}{2} e^{-2\lambda}$, and the mean of $3\lfloor(X/3)\rfloor$ equals $\lambda -1+e^{-3 \lambda /2} \left(\frac{\sin \left(\frac{\sqrt{3} \lambda }{2}\right)}{\sqrt{3}}+\cos \left(\frac{\sqrt{3} \lambda}{2}\right)\right)$:

The means for $m=1,2,3$ are shown in blue, red, and yellow, respectively, as functions of $\lambda$: asymptotically, the mean drops by $(m-1)/2$ compared to the original Poisson mean.
Similar formulas for the variances can be obtained. (They get messy as $m$ rises and so are omitted. One thing they definitively establish is that when $m \gt 1$ no multiple of $Y$ is Poisson: it does not have the characteristic equality of mean and variance)  Here is a plot of the variances as a function of $\lambda$ for $m=1,2,3$:

It is interesting that for larger values of $\lambda$ the variances increase. Intuitively, this is due to two competing phenomena: the floor function is effectively binning groups of values that originally were distinct; this must cause the variance to decrease. At the same time, as we have seen, the means are changing, too (because each bin is represented by its smallest value); this must cause a term equal to the square of the difference of means to be added back.  The increase in variance for large $\lambda$ becomes larger with larger values of $m$.
The behavior of the variance of $mY$ with $m$ is surprisingly complex. Let's end with a quick simulation (in R) showing what it can do.  The plots show the difference between the variance of $m\lfloor X/m \rfloor$ and the variance of $X$ for Poisson distributed $X$ with various values of $\lambda$ ranging from $1$ through $5000$. In all cases the plots appear to have reached their asymptotic values at the right.
set.seed(17)
par(mfrow=c(3,4))
temp <- sapply(c(1,2,5,10,20,50,100,200,500,1000,2000,5000), function(lambda) {
  x <- rpois(20000, lambda)
  v <- sapply(1:floor(lambda + 4*sqrt(lambda)), 
              function(m) var(floor(x/m)*m) - var(x))
  plot(v, type="l", xlab="", ylab="Increased variance", 
       main=toString(lambda), cex.main=.85, col="Blue", lwd=2)
})


A: As Michael Chernick says, if the individual random variables are independent then the the sum is Poisson with parameter (mean and variance) $\sum_{i=1}^{n} \lambda_i$ which you might call $\lambda$.
Dividing by $n$ reduces the mean to $\lambda / n$ and variance $\lambda / n^2$ so the variance will be less than the equivalent Poisson distribution. As Michael says, not all values will be integers.
Using the floor function reduces the mean slightly, by about $\frac12 -\frac{1}{2n}$, and affects the variance slightly too though in a more complicated manner.  Although you have integer values, the variance will still be substantially less than the mean and so you will have a narrower distribution than the Poisson. 
A: The probability mass function of the average of $n$ independent Poisson random variables can be written down explicitly, though the answer might
not help you very much.  As Michael Chernick noted in comments
on his own answer, 
the sum $\sum_i X_i$ of independent Poisson random variables $X_i$ with 
respective parameters $\lambda_i$ is a Poisson random variable with parameter $\lambda = \sum_i \lambda_i$.  Hence, 
$$P\left\{ \sum_{i=1}^n X_i= k\right\} = \exp(-\lambda)\frac{\lambda^k}{k!},
~~ k = 0, 1, 2, \ldots,$$
Thus, $\hat{Y} = n^{-1} \sum_{i=1}^n X_i$ is a random variable taking on value $k/n$ with probability $\exp(-\lambda)\frac{\lambda^k}{k!}$.  Note that
$\hat{Y}$ is not an
integer-valued random variable (though it does take on uniformly-spaced 
rational values). It follows easily that 
$Y = \lfloor \hat{Y} \rfloor$ is an integer-valued random variable
taking on value $m$ with probability
$$P\{Y = m\} = P\left\{\left\lfloor \frac{1}{n}\sum_{i=1}^n X_i
\right\rfloor = m\right\}
= \exp(-\lambda)\sum_{i=0}^{n-1}\frac{\lambda^{mn+i}}{(mn+i)!},
~~ m = 0, 1, 2, \ldots,$$
This is not the probability mass function of a Poisson random
variable.  Formulas for the mean and variance can be written down using this
probability mass function, but they don't obviously
lead to nice simple answers in terms of $\lambda$ and $n$.
Approximate values can be obtained as pointed out by Henry.
