Variance of scaled Poisson distribution

Question

I teach statistics in a community college. In my class I gave an example of image noise reduction in cosmology. For example, an object is exposed for time T and we got image $C_1(i,j)$. $C_1$ is the counts of pixel $(i,j)$ in the camera. With the lens of the camera covered, we got image $C_2(i,j)$ in a longer time of $a\times T$. Both $C_1$ and $C_2$ follow Poisson distribution. The result image after processing is $C=C_1-\frac{C_2}a$. Now my question is what the variance of $C$ should be?

I think it should be $var(C)=var(C_1)+\frac{var(C_2)}{a^2}$. But one of my students said since the count follows Poisson distribution, $var(\frac{C_2}a)$ should be equal to $\frac{C_2}a$, not $\frac{C_2}{a^2}$. It sounds weird to me, because if it is like that, $var(C)$ will be almost a constant no matter how long you take the measurements, and there will be no more benefit with long measurements. A scaled Poisson distribution is not Poisson anymore. But I cannot persuade him. Is my logic correct or what did I miss? Can anybody give some other example to explain this problem? Many thanks!

---

Thanks for your clear explanation! (@ whuber) I fully agree that the variance should be estimated by regular uncertainty propagation rules. However, I guess my student's idea is that since the process of noise generation follows Poisson distribution, intrinsically, when we get a mean of $\frac{a\nu}{a}$, we also get a variance of the same value due to the characteristics of the process. It is different from other cases such as measuring the weight of an apple. The uncertainty may depend on the smallest unit of the scale and how many times it was measured, independent of the weight of the apple. In Poisson, the mean and variance are linked and identical. But this conflicts with traditional uncertainty propagation rules and scaling of distribution. I don't know whether his view is correct or not.

Answer 1

Fix a pixel. Expose it for a unit of time, counting the incident photons. The response is the sum of two values. One is a random variable $X$ that is assumed to have a Poisson distribution. (This isn't quite true, but often is a good approximation.) Its parameter $\lambda$ is the sum of two things:

1. Something proportional to the intensity of incident radiation, multiplied by the exposure time. (This means the sensor responds linearly to the radiation, also an approximation.)

2. A term accounting for "noise," also multiplied by the exposure time. The noise is assumed to be an independent process, such as radiation originating within the instrument. (This isn't how all the noise arises, but we're ignoring that.)

Mathematically, we may express this relationship as

$$\lambda = \mu + \nu$$

where $\mu$ is the expected amount of radiation and $\nu$ (hopefully much smaller than $\mu$) is the expected amount of noise. There's no need to multiply by exposure time because it was assumed to be one unit of time.

To assess the noise, a separate exposure is made for $a$ units of time. Now $\mu$ is absent because the lens cover completely blocks the radiation (and presumably does not introduce any radiation itself!). The response $Y$ is presumed to be a Poisson variable with the same noise rate per unit of time, implying its parameter is $a\nu.$

Using properties of the Poisson distribution, we may compute the first few moments of these variables:

$$E[X] = \mu + \nu;\quad \operatorname{Var}(X) = \mu + \nu;$$

$$E[Y] = a\nu;\quad \operatorname{Var}(Y) = a\nu.$$

The expression $X-Y/a$ is an effort to adjust the image for the noise.

Linearity of expectation implies

$$E[X-Y/a] = E[X] - E[Y]/a = (\mu + \nu) - (a\nu)/a = \mu,$$

which is the signal you wish to recover: that is, $X-Y/a$ is an unbiased estimate of $\mu.$

The figure documents the results of 10,000 pixel values generated by this process with $\lambda=4,$ $\nu=2,$ and $a=10.$ The "Corrected Sample" histogram and the distribution of point heights in the "Original vs. Corrected" scatterplot clearly show that neither the corrections nor the corrected values have Poisson distributions, because they are not confined to non-negative integral values.

How far could this estimate deviate from the true value? To answer that, we need some sense of the variability of $X-Y/a.$ Its variance is easy to find and will serve well for that purpose. The independence of $X$ and $Y$ implies their covariance is zero, whence the variance of this adjusted signal is the sum of the variances of its component terms,

$$\operatorname{Var}(X - Y/a) = \operatorname{Var}(X) + \left(\frac{-1}{a}\right)^2\operatorname{Var}(Y) = (\mu + \nu) + \frac{1}{a^2}(a\nu) = \mu + \left(1 + \frac{1}{a}\right)\nu.$$

These calculations apply no matter how $X$ and $Y$ might be distributed. They rely purely on properties of variances (and some algebra) and just plug in the variances found previously. Thus, they provide useful guidance even when the many physical assumptions needed to derive them are not quite correct.

---

Notice how the variance of the noise-corrected signal $X-Y/a$ is greater than the variance of the uncorrected signal. However, when $a$ is large the $1/a$ term is relatively small, so there will be practically no difference between the two variances.

The moral of this analysis is that for the purpose of correcting a single image it's scarcely worth making the dark-lens exposure ($Y$) much longer than the original exposure. (When using $Y$ to correct many images, though, you will want to lengthen $a$ because the variance of $Y/a$ affects all the images.)