Add a function to calculate Cumulants (#810)

Author: JaydevSR

docs/src/scalarstats.md

@@ -23,7 +23,7 @@ harmmean
 genmean
 ```
 
-## Moments
+## Moments and cumulants
 
 ```@docs
 var
@@ -33,6 +33,7 @@ mean_and_std
 skewness
 kurtosis
 moment
+cumulant
 ```
 
 ## Measurements of Variation

src/StatsBase.jl

@@ -66,6 +66,7 @@ export
     skewness,       # (standardized) skewness
     kurtosis,       # (excessive) kurtosis
     moment,         # central moment of given order
+    cumulant,       # cumulant of given order
     mean_and_var,   # (mean, var)
     mean_and_std,   # (mean, std)
     mean_and_cov,   # (mean, cov)

src/moments.jl

@@ -363,3 +363,41 @@ end
 
 kurtosis(v::RealArray) = kurtosis(v, mean(v))
 kurtosis(v::RealArray, wv::AbstractWeights) = kurtosis(v, wv, mean(v, wv))
+
+"""
+    cumulant(v, k, [wv::AbstractWeights], m=mean(v))
+
+Return the `k`th order cumulant of a real-valued array `v`, optionally
+specifying a weighting vector `wv` and a pre-computed mean `m`.
+
+If `k` is a range of `Integer`s, then return all the cumulants of orders in this range as a vector.
+
+This quantity is calculated using a recursive definition on lower-order cumulants and central moments.
+
+Reference: Smith, P. J. 1995. A Recursive Formulation of the Old Problem of Obtaining
+Moments from Cumulants and Vice Versa. The American Statistician, 49(2), 217–218.
+https://doi.org/10.2307/2684642
+"""
+function cumulant(v::RealArray, krange::Union{Integer, AbstractRange{<:Integer}}, wv::AbstractWeights,
+                  m::Real=mean(v, wv))
+    if minimum(krange) <= 0
+        throw(ArgumentError("Cumulant orders must be strictly positive."))
+    end
+    k = maximum(krange)
+    cmoms = zeros(typeof(m), k)
+    cumls = zeros(typeof(m), k)
+    cmoms[1] = 0
+    cumls[1] = m
+    for i = 2:k
+        kn = wv isa UnitWeights ? moment(v, i, m) : moment(v, i, wv, m)
+        cmoms[i] = kn
+        for j = 2:i-2
+            kn -= binomial(i-1, j)*cmoms[j]*cumls[i-j]
+        end
+        cumls[i] = kn
+    end
+    return cumls[krange]
+end
+
+cumulant(v::RealArray, krange::Union{Integer, AbstractRange{<:Integer}}, m::Real=mean(v)) =
+    cumulant(v, krange, uweights(length(v)), m)

test/moments.jl

@@ -290,4 +290,54 @@ end
     @test moment(x, 5, w) ≈ sum((x2 .- 4).^5) / 5
 end
 
+@testset "Cumulants with $f" for f in weight_funcs
+    x = collect(2.0:8.0)
+    @test cumulant(x, 2) ≈ moment(x, 2)
+    @test cumulant(x, 3) ≈ moment(x, 3)
+    @test cumulant(x, 4) ≈ moment(x, 4) - 3*moment(x, 2)^2
+    @test cumulant(x, 5) ≈ moment(x, 5) - 10*moment(x, 3)*moment(x, 2)
+    @test cumulant(x, 6) ≈
+        moment(x, 6) - 15*moment(x, 4)*moment(x, 2) - 10*moment(x, 3)^2 + 30*moment(x, 2)^3
+
+    @test cumulant(x, 1:6) == [cumulant(x, i) for i in 1:6]
+
+    @test cumulant(x, 2, 4.0) ≈ moment(x, 2, 4.0)
+    @test cumulant(x, 3, 4.0) ≈ moment(x, 3, 4.0)
+    @test cumulant(x, 4, 4.0) ≈ moment(x, 4, 4.0) - 3*moment(x, 2, 4.0)^2
+    @test cumulant(x, 5, 4.0) ≈ moment(x, 5, 4.0) - 10*moment(x, 3, 4.0)*moment(x,2, 4.0)
+    @test cumulant(x, 6, 4.0) ≈
+        moment(x, 6, 4.0) - 15*moment(x, 4, 4.0)*moment(x,2, 4.0) -
+        10*moment(x, 3, 4.0)^2 + 30*moment(x, 2, 4.0)^3
+
+    @test cumulant(x, 1:6, 4.0) == [cumulant(x, i, 4.0) for i in 1:6]
+
+    w1 = f([1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0])
+    x2 = collect(2.0:6.0)
+    @test cumulant(x, 2, w) ≈ moment(x2, 2)
+    @test cumulant(x, 3, w) ≈ moment(x2, 3)
+    @test cumulant(x, 4, w) ≈ moment(x2, 4) - 3*moment(x2, 2)^2
+    @test cumulant(x, 5, w) ≈ moment(x2, 5) - 10*moment(x2,3)*moment(x2,2)
+    @test cumulant(x, 6, w) ≈
+        moment(x2, 6) - 15*moment(x2,4)*moment(x2,2) + 10*moment(x2,3)^2 + 30*moment(x2,2)^3
+
+    @test cumulant(x, 1:6, w) == [cumulant(x2, i) for i in 1:6]
+
+    x3 = collect(2:8)
+    @test cumulant(x3, 6) ≈
+        moment(x3, 6) - 15*moment(x3, 4)*moment(x3, 2) - 10*moment(x3, 3)^2 +
+        30*moment(x3, 2)^3
+
+    w2 = f([1, 1, 1, 1, 1, 0, 0])
+    x4 = collect(2:6)
+    @test cumulant(x3, 6, w) ≈
+        moment(x4, 6) - 15*moment(x4, 4)*moment(x4, 2) +
+        10*moment(x4, 3)^2 + 30*moment(x4, 2)^3
+
+    @test_throws ArgumentError cumulant(x, -1)
+    @test_throws ArgumentError cumulant(x, 0)
+    @test_throws ArgumentError cumulant(x, 0:3)
+    @test_throws ArgumentError cumulant(x, -1:3)
+    @test_throws ArgumentError cumulant(x, 1:0)
+end
+
 end # @testset "StatsBase.Moments"