# \[ANN\] BiweightStats.jl v0.2: Robust statistics based on the biweight transform

**URL:** <https://discourse.julialang.org/t/ann-biweightstats-jl-v0-2-robust-statistics-based-on-the-biweight-transform/83058>\
**Category:** Package Announcements\
**Tags:** statistics\
**Created:** [June 20, 2022, 11:23am UTC](https://discourse.julialang.org/t/ann-biweightstats-jl-v0-2-robust-statistics-based-on-the-biweight-transform/83058 "2022-06-20T11:23:18Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mileslucas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mileslucas/32/11747_2.png) [@mileslucas](https://discourse.julialang.org/u/mileslucas)\
**Post date:** [June 20, 2022, 11:23am UTC](https://discourse.julialang.org/t/ann-biweightstats-jl-v0-2-robust-statistics-based-on-the-biweight-transform/83058/1 "2022-06-20T11:23:18Z")

</div>

I’m announcing a new package, [BiweightStats.jl](https://github.com/mileslucas/BiweightStats.jl) ([docs](http://mileslucas.com/BiweightStats.jl/stable)).

It implements a suite of robust statistics based on the [biweight transform](https://www.itl.nist.gov/div898/software/dataplot/refman2/ch2/biweight.pdf). This transform, as implemented, will skip any non-finite values, and will filter any values greater than a given cutoff value from the median of the data.

The motivation for this package was for use in image processing, specifically for astronomical imaging where extreme outliers are not uncommon (e.g., from cosmic rays on CCDs). The statistics implemented are based on the NIST definitions (e.g., [scale](https://www.itl.nist.gov/div898/software/dataplot/refman2/auxillar/biwscale.htm)), and match the functionality of the [astropy.stats](https://github.com/astropy/astropy) biweight statistics.

The code has been [benchmarked](http://mileslucas.com/BiweightStats.jl/stable/#Benchmarks) against the astropy implementation, and while it readily outperforms python for most of benchmarks, I’m curious how the Julia code can be improved. As of right now, the package has 0 dependencies, but does include a hand-written iterator and loops. In theory there could be speedup by using Transducers.jl, Floops.jl, LoopVectorization.jl, Polyster.jl, etc. I believe the predominate bottlenecks are the allocations and medians required for[implementing the biweight transform](https://github.com/mileslucas/BiweightStats.jl/blob/88324b853567efe5f2a0d827b3213ad91007244a/src/BiweightStats.jl#L107-L116). I’m happy to entertain any attempts at improving the performance!

---

<div class="post-metadata">

**Author:** ![giordano](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/giordano/32/2166_2.png) [@giordano](https://discourse.julialang.org/u/giordano)\
**Post date:** [June 20, 2022, 1:18pm UTC](https://discourse.julialang.org/t/ann-biweightstats-jl-v0-2-robust-statistics-based-on-the-biweight-transform/83058/2 "2022-06-20T13:18:58Z")

</div>

I think you can save a temporary array:

```julia
julia> using Statistics, BenchmarkTools

julia> function orig(X; c=9, M=nothing)
           data = filter(isfinite, X)
           if isnothing(M)
               med = median(data)
           else
               med = M
           end
           mad = median!(abs.(data .- med))
       end
orig (generic function with 1 method)

julia> function test(X; c=9, M=nothing)
           data = filter(isfinite, X)
           med = isnothing(M) ? median(data) : M
           data .= abs.(data .- med)
           mad = median!(data)
       end
test (generic function with 1 method)

julia> @benchmark orig(X) setup=(X=randn(1000, 1000))
BenchmarkTools.Trial: 164 samples with 1 evaluation.
 Range (min … max): 21.039 ms … 35.319 ms ┊ GC (min … max): 7.46% … 0.00%
 Time (median): 27.269 ms ┊ GC (median): 1.20%
 Time (mean ± σ): 27.172 ms ± 2.476 ms ┊ GC (mean ± σ): 2.73% ± 3.26%

                    █▃▁▄▁ ▃▄▁▃ ▁▃▄▃▁▁▃▃▆ ▆ ▃ ▁       
  ▄▁▁▄▁▁▁▁▁▆▄▇▇▆▇▆▄▆█████▇▆████▇▆█████████▇▇█▇▆▄█▆▇▇▁▄▁█▄▁▁▁▄ ▄
  21 ms Histogram: frequency by time 32.5 ms <

 Memory estimate: 22.89 MiB, allocs estimate: 6.

julia> @benchmark test(X) setup=(X=randn(1000, 1000))
BenchmarkTools.Trial: 183 samples with 1 evaluation.
 Range (min … max): 18.818 ms … 31.551 ms ┊ GC (min … max): 0.00% … 0.95%
 Time (median): 25.092 ms ┊ GC (median): 0.00%
 Time (mean ± σ): 24.793 ms ± 2.497 ms ┊ GC (mean ± σ): 1.28% ± 2.01%

               ▃ ▂ ▅ ▅ ▂ ▂▂ ▅█▆▃▅▂▃ ▃▃ ▂ ▂          
  ▄▁▄▁▄▄▄▁▅▇▄▄▇█▅▇▇███████▅████▅▅██████████▅█▇▄▅▄▅▄▄█▁▄▄▁▁▄▁▅ ▄
  18.8 ms Histogram: frequency by time 30.6 ms <

 Memory estimate: 15.26 MiB, allocs estimate: 4.

```

It’s a ~10% improvement in this run.

---

<div class="post-metadata">

**Author:** ![mileslucas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mileslucas/32/11747_2.png) [@mileslucas](https://discourse.julialang.org/u/mileslucas)\
**Post date:** [June 20, 2022, 6:37pm UTC](https://discourse.julialang.org/t/ann-biweightstats-jl-v0-2-robust-statistics-based-on-the-biweight-transform/83058/3 "2022-06-20T18:37:58Z")

</div>

Thanks @giordano! I played around with this a little before realizing it caused a bug with the covariance!

We can’t reorder the elements when calculating the covariance between two vectors. Well, they can be reordered, just so long as _both_ vectors are reordered the same way. This means I can’t use `median!` on the intermediate data for that method.

What I can do is add a keyword argument whether its safe to do this, but it means the `midcov` and `midcor` methods won’t see the speedup.