examples (c, local n most recent observations) from the redesign of RollingFunctions.jl

ShowAndTell(c).jl

# see example at line 80

using StatsBase

function StatsBase.median(a::T, b::T, c::T) where {T}
    b, c = minmax(b, c)
    a, c = minmax(a, c)
    max(a, b)
end

abstract type Accumulator{T} <: Function end

const AccNum = Float64 # client settable via ENV

StatsBase.nobs(acc::Accumulator) = acc.nobs

fn(acc::Accumulator) = acc.fn

mem(acc::Accumulator) = acc.mem
mem(acc::Accumulator, idx::Int) = @inbounds mem(acc)[idx]
mem(acc::Accumulator, idxs::StepRange) = @inbounds mem(acc)[idxs]

nmem(acc::Accumulator) = length(mem(acc))


#=
   incremental accumulation of the mean and standard deviation
   using the median of the current and two most recent observations 
=#

mutable struct AccMeanStd{T,N} <: Accumulator{T}
    nobs::Int
    mean::T
    svar::T
    mem::NTuple{N,T}
end

function AccMeanStd(::Type{T}=AccNum, nmem::Int=2) where {T}
     AccMeanStd{T,nmem}(0, zero(T), zero(T), ntuple(z->zero(T), nmem))
end

function (acc::AccMeanStd{T})() where {T}
    if isone(acc.nobs)
       (acc.mem[1], zero(T))
    else
        unbiased_var = acc.svar / (acc.nobs- 1)
        (acc.mean, sqrt(unbiased_var))
    end    
end

function (acc::AccMeanStd{T})(x) where {T}
    acc.nobs += 1    
    prior_mean = acc.mean

    y = median(x, acc.mem[1], acc.mem[2])
    change = y - prior_mean

    acc.mean = change / acc.nobs + prior_mean
    acc.svar = change * (y - acc.mean) + acc.svar

    acc.mem = (x, acc.mem[1])
    acc
end

function (acc::AccMeanStd{T})(xs::A) where {T, A<:AbstractVector{T}}
    for i in eachindex(xs)
        acc(xs[i])
    end
    acc
end

function (acc::AccMeanStd{T})(xs::NTuple{N,T}) where {T, N}
    for i in eachindex(xs)
        acc(xs[i])
    end
    acc
end






#= 
   the following example shows
 
  mean of median of 3 at index 1 = 5.0
 standard deviation  at index 1 = 0.0

 mean of median of 3 at index 2 = 2.5
 standard deviation  at index 2 = 3.535534

 mean of median of 3 at index 3 = 3.3333333
 standard deviation  at index 3 = 2.8867514

 mean of median of 3 at index 4 = 3.5
 standard deviation  at index 4 = 2.3804762

 mean of median of 3 at index 5 = 2.2
 standard deviation  at index 5 = 3.5637062

 mean of median of 3 at index 6 = 1.3333335
 standard deviation  at index 6 = 3.8297086

 mean of median of 3 at index 7 = 1.2857144
 standard deviation  at index 7 = 3.498299

 mean of median of 3 at index 8 = 1.2500001
 standard deviation  at index 8 = 3.2403703

=#

xs = Float32[5, 8, 4, -3, -6, 3, 1, -9]

registers = 2
acc = AccMeanStd(eltype(xs), registers)

for i in eachindex(xs)
    acc(xs[i])
    mean_std = acc()
    
    str = string(" mean of median of 3 at index ", i, " = ", mean_std[1])
    println(str)
    str = string(" standard deviation  at index ", i, " = ", mean_std[2])
    println(str)
    println("")
end