# Construct probability that measurement belongs to one of two sets

We are dealing with a measurement apparatus that is cursed by noise and are trying to find out if a measurement was noise or an actual measurement.

Assume we have a beam of light incident on an square array of photo detectors. The "counts" measured in the individual tubes follow Poisson statistics. A typical beam distributes light across 4 photo tubes, most in the tube closest to the point of incident. Noise would be distributed uncorrelated.

We can measure the count distributions for a location where we know we had an actual beam hitting, and for a location where we know no beam was present. From these measurements we can construct empirical count PDFs for measurements in the individual tubes, $p_{i,\,\rm{true}}$ and $p_{i,\,\rm{noise}}$.

Since for a noise measurement the counts in the individual tubes are uncorrelated I would guess that

$$P_\mathrm{noise} = \prod_i p_{i,\,\mathrm{noise}}$$

would be a good guess for the probability that this measurement was noise. Can I use a similar expression

$$P_\mathrm{true} = \prod_i p_{i,\,\mathrm{true}}$$

even though counts would be correlated for the "true" case?

• What is a "measurement"? Is it (a) a determination that the array has detected a signal; (b) a determination for each pixel in the array as to whether or not it has detected a signal; (c) an estimate of total number of photons reaching the array; (d) an estimate of photons reaching each pixel (i.e., an image); or (e) something else?
– whuber
Feb 22, 2011 at 18:40
• @whuber: We are primarily interested to know (a) if we measured a signal, but at a later stage would be interested in (c) as well. Feb 22, 2011 at 18:43
• Is this a series of discrete measurements or is the array constantly being monitored? In other words, can we treat this as a hypothesis test or is it a control chart/changepoint/signal detection problem?
– whuber
Feb 22, 2011 at 20:11
• @whuber: We take measurements at singular points in time. Feb 22, 2011 at 20:16

## 1 Answer

According to the comments after the question, this is a hypothesis testing situation. You have stipulated that you can accurately assess the null distribution of the individual cell counts. We need a test statistic. The nature of the problem suggests running a small kernel over the array (essentially to deconvolve the signal). A simple choice would be a 2 x 2 mean. We would then look at the maximum of those values. The Neyman-Pearson lemma says that a good test is based on setting a critical value for this statistic. The test size depends on the null distribution of the test statistic while its power depends on the intensity of the signal. You need to choose a critical value to appropriately balance the expected false positive rate with the power for the kinds of signals you are looking for.

All this is routine; the only non-routine aspect is determining the null distribution of the statistic. It's the maximum of a set of fairly highly correlated linear combinations of Poisson variates. Specifically, if we index the rows and columns of the array and let $X_{i,j}$ be the Poisson variate for the cell in row $i$, $1 \le i \le m$, and column $j$, $1 \le j \le n$, then the 2 x 2 mean is the array of values $Y_{i,j} = (X_{i,j} + X_{i+1,j} + X_{i,j+1} + X_{i+1,j+1})/4$, $1 \le i \le m-1$, $1 \le j \le n-1$, and the test statistic is $t = \max\{Y_{i,j}\}$.

If someone could jump in and tell us how to compute the distribution of $t$ that would be nice, but I suspect it's not an easy calculation. Would you consider a small simulation? It's easy to set up in R or Mathematica, for instance.

Alternatively, use an approximation. The results of a few simulations with arrays from 9 to 1600 elements and photon intensities from 1 to 4096 per cell are consistent with two obvious approximations: one can treat all $(m-1)(n-1)$ of the $Y_{i,j}$ as independent or one can take every fourth one and treat them as independent, and then calculate the distribution of their maximum. For small intensities the upper tails of the simulated distributions of $t$ (10,000 iterations per simulation) appear to behave like the latter approximation: that is, $t$ behaves like the largest of $(m-1)(n-1)/4$ independent averages of four Poisson variates. For large intensities the tails are quite close to the former approximation: that is, $t$ behaves like the largest of $(m-1)(n-1)$ independent averages of Poisson variates.

This is a histogram of $t$ for a 30 by 40 grid with cell intensity 100 (10,000 simulations). Because its 99th quantile equals 122.25, you could create a test with at most 1% false positive rate by setting the critical value to 123. An average photon count of around 135 in a single block of four cells would be readily detectable with this method. That would represent a total flux of 4*(135 - 100) = 140 photons above background.

For the record, here is the Mathematica code used to generate this histogram.

simulate[m_Integer, n_Integer, \[Mu]_] /;
m >= 2 && n >= 2 && \[Mu] > 0 := With[
{f = PoissonDistribution[\[Mu]]},
y = ListConvolve[{{1, 1}, {1, 1}}/4, RandomInteger[f, {m, n}]]
];
With[{m = 30, n = 40, \[Mu] = 100, nTrials = 10000},
null = ParallelTable[Max[Flatten[simulate[m, n, \[Mu]]]], {i, 1, nTrials}];
Histogram[null, {1/4}, AxesLabel -> {"t", "Count"}]
]

• @whuber (+1) I think there's a typo with the closing ] around function args list (should be before :=).
– chl
Feb 23, 2011 at 15:39
• @chl You have great eyes! However, this was cut and pasted from working code. Mathematica is flexible about where the criteria (following "/;") need to appear. In my experience they are most reliably interpreted when they all appear outside the arguments list, as shown here, even though in some cases that can slow down execution slightly.
– whuber
Feb 23, 2011 at 15:46
• @whuber Indeed, it works if there's no carriage return before /; -- my fault!
– chl
Feb 23, 2011 at 15:55
• @chl No, you're right. I had changed the carriage returns to fit the pasted version within the space allotted. It turns out Mathematica objects when a carriage return precedes /; but not if it follows /;. (I call that bad parsing, but that's another matter.) I have changed it so that anyone who pastes this code wholesale back into Mathematica will encounter no error messages.
– whuber
Feb 23, 2011 at 16:01
• Thank you for this detailed answer whuber. You didn't however not address my question (last sentence) at all. Could you make a small comment (in your answer) on that, please? Mar 1, 2011 at 17:30