Smoothing Noisy Data using Moving Mean

General Usage
,
1

Hi there,

Is there a function in Julia that is similar to MATLAB’s smooth() function.

MATLAB smooth function: z = smooth(y, span, method) where

y → input array;
span → span of moving average;
method → option between default moving average, lowess, loess etc.

THe result of this operation z is a vector of the same dimension as y which is basically a smoothening of the original array y .

2 Likes
2

have you seen RollingFunctions.jl? Seems to be close to what you want. But it doesn’t have out-of-the-box support for loess etc. It does have a way to use a custom function, so maybe it’s easy to add on your own.

1 Like
3

For LOESS,

https://github.com/JuliaStats/Loess.jl

1 Like
4

If it’s just a moving average and you want the same exact results as MATLAB then the only thing to figure out is what strategy MATLAB uses to estimate the tails (AR models, … etc.)

If you just want a smoother then they all have pros and cons, e.g. some can interpolate some cannot.

Anyway, to add options to the ones above we also have a Henderson moving average filter here GitHub - viraltux/Forecast.jl: Julia package containing utilities intended for Time Series analysis.

using Plots, Forecast
x = log.( collect(1:.1:11) .+ rand(101))

scatter(x)
plot!(hma(x,21),legend=nothing)

image

2 Likes
5

There’s no singular function (yet) that has parity with MATLAB’s smoothdata(), and I +1 for whoever would like to implement something like it!
This works for a moving average at least:

using NaNStatistics, Plots

y = randn(1000)
window_size = length(y)/100.0
z = movmean(y, window_size)
plot(y, linestyle = :dot)
plot!(z, linewidth = 2)

smooth

2 Likes
6

If you just want to plot the data more smoothly, there is this one:

julia> using EasyFit

julia> x = rand(100);

julia> m10 = movavg(x,10)

 ------------------- Moving Average ---------- 

 Number of points averaged: 11 (± 5 points)

 Pearson correlation coefficient, R = 0.43442913764734875

 Averaged X: x = [0.4737608452310648, 0.5121944744036446...
 residues = [-0.3269117547107911, 0.40109456416787137...

 -------------------------------------------- 

julia> m50 = movavg(x,50)

 ------------------- Moving Average ---------- 

 Number of points averaged: 51 (± 25 points)

 Pearson correlation coefficient, R = -0.014175383510694936

 Averaged X: x = [0.5225594827160708, 0.5220326497692106...
 residues = [-0.37571039219579705, 0.3912563888023054...

 -------------------------------------------- 

julia> plot([x,m10.x,m50.x],linewidth=2)

image

7

For a package-free moving average, see @tim.holy’s solution here, adapted as follows for m odd integer > 1:

function moving_average(A::AbstractArray, m::Int)
    out = similar(A)
    R = CartesianIndices(A)
    Ifirst, Ilast = first(R), last(R)
    I1 = m÷2 * oneunit(Ifirst)
    for I in R
        n, s = 0, zero(eltype(out))
        for J in max(Ifirst, I-I1):min(Ilast, I+I1)
            s += A[J]
            n += 1
        end
        out[I] = s/n
    end
    return out
end
4 Likes
8

ImageFiltering.jl has multidmensional array smoothing. There’s also a mapwindow(f, A, window) which allows you to apply f to “windows” of A rather than the single elements that map(f, A) applies f to.

2 Likes
9

I’ll also note that an exponentially weighted mean is similar to a moving window and is easy to use via OnlineStats:

o = Mean(weight = ExponentialWeight(.1))

exp_smooth = [value(fit!(o, yi)) for yi in y]
1 Like
10

Hi,
How do you choose the optimal window size?
The package is very promising, I hope it develops quickly, I am looking forward to using it soon.
Many thanks,
Ayush