[Python] Why sigma clipping need iterations?

If I only want the high-SNR data, I do sigma-clipping to an array.

As this link says

Suppose you have a set of data. Compute its median m and its standard deviation sigma. Keep only the data that falls in the range (m-a*sigma,m+a*sigma) for some value of a, and discard everything else.

Then it just need one step to get the (m-a*sigma,m+a*sigma) part:

lower_bound = m - sigma * a
upper_bound = m + sigma * a

mask = (data >= lower_bound) & (data <= upper_bound)
cleaned_data = data[mask]


Why there is a need for iterating?


Then I do a test upon a very large (1000 * 800 * 600) astronomy data cube

from astropy.stats import sigma_clip
i = 1
clipped_data = sigma_clip(data, sigma=3, maxiters=i, masked=True)
mask = clipped_data.mask # data outside 3σ
N = np.count_nonzero(mask)
print(N)


then I got this:

i = 1 N = 4854040
2 6930075
3 7541519
4 7715989
5 7765465
7 7783735
9 7785237




On the other hand, if I do this in pure numpy:

m = np.median(data)
sigma = np.std(data)
a = 3
lower_bound = m - sigma * a
upper_bound = m + sigma * a

mask = (data <= lower_bound)
mask2 = (data >= upper_bound)
np.count_nonzero(mask+mask2)


This gives 4854042, close to the result given by iteration = 1

So why when I add iterations, and the number of data points outside 3σ will be more?

Continue reading...