These are chat archives for astropy/astropy

8th
Nov 2018
Till Stensitzki
@Tillsten
Nov 08 2018 22:17
Hi
I may be wrong, but I think the current sigma_clip function is buggy.
I wrote a simple version of the clipping function myself:
Till Stensitzki
@Tillsten
Nov 08 2018 22:23
def sigma_clip(data, sigma=3, iterations=5, axis=-1):
    data = np.ma.masked_invalid(data)
    num_masked = 0
    for i in range(iterations):
        median = np.ma.medan(data, axis, keepdims=1)
        std = np.std(data, axis, keepdims=1)
        upper, lower = median + sigma*std, median - sigma*std
        data = np.ma.masked_greater(data, upper, copy=False)
        data = np.ma.masked_less(data, lower, copy=False)
        n = data.mask.sum()
        if n == num_masked:
            break
        else:
            num_masked = n
    return data
When I now compare both functions the results are slightly different, but astopys sigma_clip seems to mask a little too much.
This is visible in the following example:
np.random.seed(0)
test_data = np.random.randn(100)
test_data[np.random.randint(0, 100, 10)] *= 10
plt.plot(test_data, '.')
plt.plot(sigma_clip(test_data, 2), 'o')
plt.plot(st.sigma_clip(test_data, 2), 'o')
Till Stensitzki
@Tillsten
Nov 08 2018 22:36
Oh not is not, sorry for that. The bug is not visible for one-dimensional data
Till Stensitzki
@Tillsten
Nov 08 2018 22:47
Ok, luckily I am the slow one here, I taught astropy axis defaults to -1. Instead None defaults to clipping the flat array. Sorry for the noise.