Guide

Choosing an output

The thresholds are positional and have no default. Every function takes the lower and upper thresholds as sigma_lower, sigma_upper, as in sigma_clip(x, 2, 4). Passing a single sigma is shorthand for passing it twice: sigma_clip(x, 3) is the same call as sigma_clip(x, 3, 3), and the same holds for every function in the table. The thresholds are in units of the spread and must be finite and positive.

All of them accept the same keywords:

The input must be an AbstractArray. Iterators such as generators, Iterators.filter or skipmissing are not accepted; collect them first:

julia> x = [1.0, missing, 2.0, 3.0, 4.0, 100.0];

julia> sigma_clip_bounds(collect(skipmissing(x)), 3)
(-1.9478066555168052, 6.947806655516805)

Clipping in place

sigma_clip! writes NaN, so it needs a floating-point array. Here the lower threshold is 2 and the upper threshold is 4:

julia> x = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, -10.0, 12.0];

julia> sigma_clip!(x, 2, 4);

julia> x'
1×9 adjoint(::Vector{Float64}) with eltype Float64:
 0.0  1.0  2.0  3.0  4.0  5.0  6.0  NaN  12.0

For integer arrays, use sigma_clip, which returns a floating-point copy.

Excluding values from bound estimation

exclude keeps values out of the center and spread estimates. The excluded values are still compared with the final bounds, so it does not protect them from being classified as outliers:

julia> data = [-100.0, 0.0, 0.1, -0.1, 50.0, 0.05];

julia> exclude = Bool[true, false, false, false, true, false];

julia> sigma_clip_mask(data, 3; exclude)'
1×6 adjoint(::BitVector) with eltype Bool:
 0  1  1  1  0  1

Statistics on the retained values

sigma_clipped_stats always returns the final center and spread. Extra statistics are requested with symbol-keyed pairs; each function receives a vector of the retained values:

julia> using Statistics

julia> data = [1.0, 2.0, 3.0, 4.0, 5.0, 100.0];

julia> sigma_clipped_stats(data, 3; :mean => mean, :n => length, :max => maximum)
(center = 3.0, spread = 1.4826022185056018, mean = 3.0, n = 5, max = 5.0)

Write pair keys as :name => f. Without the colon, name is looked up as a variable in the calling scope. Every keyword other than workspace, exclude, center, spread and maxiter is taken as a statistic, so a misspelled keyword fails when its value is called on the retained values.

Custom center and spread

Any callable that takes an AbstractVector and returns a scalar can be used as center or spread. It may reorder its argument, which is SigmaClip's internal buffer and not the user's array:

julia> using Statistics

julia> data = [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 100.0];

julia> iqr_spread(v) = (quantile(v, 0.75) - quantile(v, 0.25)) / 1.349;

julia> sigma_clip(data, 3; spread = iqr_spread)'
1×9 adjoint(::Vector{Float64}) with eltype Float64:
 1.0  2.0  3.0  4.0  5.0  6.0  7.0  8.0  NaN

fast_median! and mad_std! can also be called directly. Both reorder their input.

Reusing a workspace

Every call needs scratch space for a copy of the finite values. When clipping many arrays, allocate it once with SigmaClipWorkspace and pass it through the workspace keyword:

julia> image = repeat([0.0 1.0 2.0 3.0 4.0 5.0 6.0 50.0], 3);

julia> n = size(image, 2);

julia> workspace = SigmaClipWorkspace(Vector{Float64}(undef, n), Vector{Float64}(undef, n));

julia> for row in eachrow(image)
           sigma_clip!(row, 3; workspace)
       end

julia> image
3×8 Matrix{Float64}:
 0.0  1.0  2.0  3.0  4.0  5.0  6.0  NaN
 0.0  1.0  2.0  3.0  4.0  5.0  6.0  NaN
 0.0  1.0  2.0  3.0  4.0  5.0  6.0  NaN

The main buffer must have the input's element type and at least as many elements as the input. The auxiliary buffer is used only by mad_std!; pass nothing when the spread estimator does not need it:

julia> using Statistics

julia> data = [1.0, 2.0, 3.0, 4.0, 5.0, 1.0, 2.0, 3.0, 4.0, 5.0, 100.0];

julia> workspace = SigmaClipWorkspace(similar(data), nothing);

julia> sigma_clip_bounds(data, 3; workspace, spread = std)
(-1.4721359549995796, 7.47213595499958)

Custom workspace types

Any type can serve as a workspace by implementing SigmaClip.workspace_buffer and SigmaClip.workspace_auxbuffer, with the same requirements on the returned buffers:

julia> struct MyWorkspace{B, A}
           buf::B
           aux::A
       end

julia> SigmaClip.workspace_buffer(ws::MyWorkspace) = ws.buf;

julia> SigmaClip.workspace_auxbuffer(ws::MyWorkspace) = ws.aux;

julia> workspace = MyWorkspace(Vector{Float64}(undef, 16), Vector{Float64}(undef, 16));

julia> sigma_clip_bounds([1.0, 2.0, 3.0, 4.0, 5.0, 100.0], 3; workspace)
(-1.4478066555168052, 7.447806655516805)