What is the variance of the least of several series? Given a number of series (or images) that are independent and have a Poisson distribution with the same mean, what is the variance of the series generated from taking the point-by-point minimum?
i.e. if the original data are I(i,j), then the resultant series is
M(i) = min(I(i,*)).

Clearly it is less than the variance of the original series, and by playing with a random number generator I estimate the factor to be about 2/3 for 2 series and 4/9 for 4, but I would like to have a proper expression for an arbitrary number of series.
 A: I will avoid the image aspect, and just concentrate on one pixel. Then we have $X_1, \dotsc, X_n$ iid variables $\mathcal{P}(\lambda)$ and interest is on $M=\min_{i=1\dotsc,n} X_i$. First, what is the distribution? Using standard order statistics arguments, we have 
$$
   P(M \le m)=1-P(X_1 > m)^n
$$ and 
$$
P(M=m) = P(X_1 > m-1)^n - P(X_1 > m)^n
$$
Exact expression for expectation and variance will not be available, maybe some analytic approximations can be found, but here I will show how to compute numerically in R.  First, some code:
dM <- function(x, lambda, n, log=FALSE) {
    upperm1.log <- n*ppois(x-1, lambda, lower.tail=FALSE, log.p=TRUE)
    upper.log   <-  n*ppois(x, lambda, lower.tail=FALSE, log.p=TRUE)
    if(log) return(log(exp(upperm1.log)-exp(upper.log))) else return(exp(upperm1.log)-exp(upper.log))
    }
pM <- function(q, lambda, n, lower.tail=TRUE, log.p=FALSE) {
    upper.log <- n*ppois(q, lambda, lower.tail=FALSE, log.p=TRUE)
    if(lower.tail) { if(log.p) return(log(1-exp(upper.log))) else return(1-exp(upper.log))} else        {if(log.p) return(upper.log) else return(exp(upper.log))}
}

Using this we can approximate expectation and variance numerically, by replacing the infinite sums by a finite sum. Some simple code:
E_M <- Vectorize(function(lambda, n) {
    upperlim <- ceiling(20*lambda) #should be enough,  but not investigated ...
    sum( pM(0:upperlim, lambda, n, lower.tail=FALSE) )
    })   

V_M <- function(lambda, n) {
    upperlim <- max(ceiling(100*lambda))
    e <- E_M(lambda, n)
    N <- length(lambda)
    result <- numeric(N)  
    for (i in seq(along=lambda)) result[i] <- sum( ((0:upperlim)-e[i])^2*dM(0:upperlim, lambda, n) )
    result  
}

An example of use:
 V_M(1:10, 5)
 [1] 30.19555 25.63541 21.16186 18.36645 17.88784 20.09637 25.24123 33.50449
 [9] 45.02656 59.91987
 V_M(1:10, 10)
 [1]  6.976930  6.362066  5.601977  5.603438  6.759195  9.305575 13.407480
 [8] 19.188868 26.747453 36.162877

It makes sense that the minimum is less variable when taken over more observations. 
