I have some more thoughts on how to make Windsorized sigma clipping work better but don't have time to get into it right now. Hopefully, I can pull some thoughts together this evening.

Tim

OK. I finally have a chance to get back to this. To put things bluntly, I think we can do pixel rejection better. I'll get to why in a minute. First a little background. Part of what I do for a living is radiation detection and counting. There is a whole lot of crossover between that and astrophotography. Both deal with Poisson noise. Both are trying to separate a desired signal from that noise. Both deal with confounding signals from other sources like system electronics, etc. So, I may have some insight from a perspective not frequently seen.

Let's get back to sigma clipping rejection. Whether Winsorized (my apologies to Mr. Winsor for misspelling his name earlier) or not, the idea is to establish a central value that gets as close as possible to the "true" mean of the actual sky brightness for that pixel, establish some clipping criteria high and low and then reject pixels outside of those criteria when averaging the stack. The clipping criteria are defined as a multiple of the standard deviation (sigma) of the stack. Both sigma clipping and Winsorized sigma clipping are iterative processes that are intended to converge on the "true" signal distribution and identify and reject outliers.

With sigma clipping (note that this is for PI algorithms only) we start out with the full pixel stack, complete with outliers high and low that are not representative of the Poisson distribution of the signal. The routine selects the median (not the mean) from the raw stack, and calculated the standard deviation of the raw stack. Right away, there is a problem that can best be illustrated by an example. Lets pick an arbitrary dataset:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000, 1500, 10000

It should be clear that the actual target distribution is centered around 500, with some clear outliers on each end. We'll run through a couple of iterations of sigma clipping on this dataset. To start, we need the median and standard deviation. The median is just the central value in the ordered set. Conveniently in this dataset, the median is 500. The standard deviation is a simple calculation and comes out to be 2156. Right away, we can see that there is a problem. The actual target distribution clearly starts at about 494 and ends up at about 506. If we do a quick standard deviation for what is clear to us is the standard deviation, we get about 4. That's a whole lot less than 2156. Clearly, the first shot at the standard deviation for sigma clipping is a terrible estimate. But let's keep going. The routine will take the median and standard deviation it just calculated (500 and 2156) and use the sigma high and low rejection numbers the users enters and goes through each pixel in the stack and tests it. We'll use 3 and 3 for the sigma high and low for this example and just test the pixels on each end, starting with 19. The test is (median - pixelvalue)/sigma for low rejection and (pixelvalue - median)/sigma for high rejection. So for the first pixel with a pixel value of 19, the calculation is (500 - 19)/2156 = 0.22. If this number is greater than 3 (sigma low) it gets rejected. Pixel one does not get rejected. Let's move on to pixel 20 with a pixel value of 10,000. The calculation is (10,000 - 500)/2156 = 4.4. This is greater than our sigma high value of 3, so it gets rejected as it clearly should. I've already done the math and no other pixels get rejected. So the new stack is:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000, 1500

The routine takes this stack and starts a new iteration. Now, the median stays the same at 500, but the new standard deviation is 346 - clearly an improvement, but still nowhere near close to 4. If we run the rejection tests again, the only pixel that gets rejected is pixel #19 with a value of 1500. The new stack is:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000

Here again, the median and sigma are calculated: median = 499, sigma = 249. Still far from the 4 or so it should be but still better than previous iterations. But, here things get interesting. When the sigma high and sigma low tests are applied, **no pixels are rejected**. So, if we're not paying very close attention, the routine will finish here. It will have rejected the two brightest pixels, so may very well have eliminated the satellite trail or hot pixel in the stack, but the 18-pixel stack still contains 5 pixels that clearly should not be included. As a result, the stack that ends up in the final image has a sigma of 249, when it should have a sigma of around 4. That's terrible. The result will be what looks like shot noise, but is in fact poor pixel rejection. Of course, we can tweak the sigma high and low until we get a picture that looks better, but that is very subjective.

Let's move on to how Winsorized sigma clipping works. Let's start with the same dataset:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000, 1500, 10000

Like before, a median and standard deviation are calculated: median = 500, sigma = 2156. Here, however, two different threshold values are calculated. The low threshold is:

Tlow = median - 1.5* sigma = -2734

And the high threshold is:

Thigh = median + 1.5*sigma = 3734

The 1.5 is a somewhat arbitrary factor hardwired into the PI routine. The routine uses these threshold values to test each pixel. If it is lower than the low threshold or higher than the high threshold, it gets **replaced** by the nearest value. So, after the first iteration, pixel #20 with a value of 10,000 get replaced with pixel #19's value of 1500. So the new stack is:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000, 1500, 1500

Then the next iteration begins with a median of 500 and a sigma of 407. Here, PI does another operation by multiplying the sigma by 1.134 - another hardwired number. The resulting sigma is 462, and Tlow and Thigh are -193 and 1193. The next stack will look like:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000, 1000, 1000

The next iteration results in no change. It is important to note here that we're not choosing the final stack with this process so far, we're just Winsorizing the set to calculate a better sigma. So, after Winsorizing the dataset, we get a median of 500 and a new sigma of 296. Note that the Winsorized sigma is actually worse than the sigma that the sigma clipping process settled on. This is an indication that 20 may be two few subs for Winsorizing to have an advantage. Clearly it does not in this case. At any rate, we now apply the sigma high and sigma low thresholds defined by the user (we used 3 and 3 above) to the **original** dataset. Original:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000, 1500, 10000

And the reduced dataset after sigma clipping:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000

Note that the final stack is identical to the sigma clipping stack. So, for this dataset, there is no difference whatsoever between sigma clipping and Winsorized sigma clipping. Furthermore, both algorithms left a large number of outliers in the stack and ended up with a stack having a very high sigma - nearly two orders of magnitude higher than it should be.

So, what can we do? Of course, we can tweak settings like I mentioned above. Tighter sigma values would certainly help, but that is a bit arbitrary. So, how does this relate to radiation detection? Well, we have a similar problem in that field. We need to separate background radiation from source radiation. The analogy to pixel rejection is that the source radiation is the outlier, and background radiation is the Poisson distribution we are trying to separate it from. Rather than applying sigma clipping rejection to identify the outliers, we use confidence intervals. The background distribution is a true Poisson distribution and we can establish confidence intervals that give us the probability of a data point being within the distribution or not. We also take advantage of the convenient fact that the standard deviation of a Poisson distribution is simply the square root of the mean. There is a relatively famous equation by Currie and Brodsky that allows us to establish whatever confidence intervals we want and all we need to know is the mean of the background and the count time - which can be represented by number of subs if the subs are the same length.

In terms of stacking subs, we can't count the background with full confidence of excluding the source like we can in radiation detection. But just as PI does, we can fairly accurately guess at the mean by selecting the median. The more subs we get, the closer the median will be to the mean, and the median is largely immune to outliers if the number of outliers is small compared to the number of data points in the target distribution. So, I'll use the median as a surrogate for the mean. The only other thing we need is the z-score for the confidence interval we desire. I think I would like to have 95% confidence that the pixels I keep are within the target distribution. The z-score for a two-tailed 95% confidence interval is 1.96. Adapting the Currie equation to our needs, I get:

Delta = 1.96*sqrt(2*N*median)/N

As you can see, delta depends only on the median and the number of subs (N). Pixel rejection would simply be those pixels that lie outside of median +/- delta. Let's see how that works for our dataset:

median = 500

N = 20

delta = 14

So our rejection thresholds are 486 low and 514 high. Let's see what that does to our dataset. Here is the original:

19, 20, 21, 22, 494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506, 1000, 1500, 10000

And here is the "confidence clipped" (I think I just coined a phrase too) dataset:

494, 495, 496, 497, 498, 500, 500 501, 501, 501, 504, 505, 506

As you can see, this simple, non-iterative approach very effectively eliminated all outliers and retained the target distribution set perfectly. We will **only** be stacking target distribution pixels and look at the sigma of our new stack - it's **3.85**. Compare that to the sigma for sigma clipping (249) and Winsorized sigma clipping (296).

As stacks get larger, the statistics only improve as you can see by the relationship that N has to the equation.

I'm going to go ahead and be a bit more blunt. I think we're doing pixel rejection **wrong.** I think this is the way to do it.

Thoughts?

Tim

**Edited by spokeshave, 06 October 2016 - 02:36 PM.**