Preserve eltype where possible for moments

src/moments.jl

@@ -42,7 +42,7 @@ function var(v::RealArray, w::AbstractWeights; mean=nothing,
                   corrected::DepBool=nothing)
     corrected = depcheck(:var, corrected)
 
-    if mean == nothing
+    if mean === nothing
         varm(v, w, Statistics.mean(v, w); corrected=corrected)
     else
         varm(v, w, mean; corrected=corrected)
@@ -85,14 +85,14 @@ end
 function varm(A::RealArray, w::AbstractWeights, M::RealArray, dim::Int;
                    corrected::DepBool=nothing)
     corrected = depcheck(:varm, corrected)
-    varm!(similar(A, Float64, Base.reduced_indices(axes(A), dim)), A, w, M,
+    varm!(similar(A, promote_type(eltype(A), eltype(w)), Base.reduced_indices(axes(A), dim)), A, w, M,
                dim; corrected=corrected)
 end
 
 function var(A::RealArray, w::AbstractWeights, dim::Int; mean=nothing,
                   corrected::DepBool=nothing)
     corrected = depcheck(:var, corrected)
-    var!(similar(A, Float64, Base.reduced_indices(axes(A), dim)), A, w, dim;
+    var!(similar(A, promote_type(eltype(A), eltype(w)), Base.reduced_indices(axes(A), dim)), A, w, dim;
          mean=mean, corrected=corrected)
 end
 
@@ -223,7 +223,7 @@ end
 ##### General central moment
 function _moment2(v::RealArray, m::Real; corrected=false)
     n = length(v)
-    s = 0.0
+    s = zero(eltype(v))
     for i = 1:n
         @inbounds z = v[i] - m
         s += z * z
@@ -233,7 +233,7 @@ end
 
 function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false)
     n = length(v)
-    s = 0.0
+    s = zero(promote_type(eltype(v), eltype(wv)))
     for i = 1:n
         @inbounds z = v[i] - m
         @inbounds s += (z * z) * wv[i]
@@ -244,7 +244,7 @@ end
 
 function _moment3(v::RealArray, m::Real)
     n = length(v)
-    s = 0.0
+    s = zero(eltype(v))
     for i = 1:n
         @inbounds z = v[i] - m
         s += z * z * z
@@ -254,7 +254,7 @@ end
 
 function _moment3(v::RealArray, wv::AbstractWeights, m::Real)
     n = length(v)
-    s = 0.0
+    s = zero(promote_type(eltype(v), eltype(wv)))
     for i = 1:n
         @inbounds z = v[i] - m
         @inbounds s += (z * z * z) * wv[i]
@@ -264,17 +264,17 @@ end
 
 function _moment4(v::RealArray, m::Real)
     n = length(v)
-    s = 0.0
+    s = zero(eltype(v))
     for i = 1:n
         @inbounds z = v[i] - m
         s += abs2(z * z)
     end
     s / n
 end
 
-function _moment4(v::RealArray, wv::AbstractWeights, m::Real)
+function _moment4(v::RealArray, wv::AbstractWeights, m::Real; corrected=false)
     n = length(v)
-    s = 0.0
+    s = zero(promote_type(eltype(v), eltype(wv)))
     for i = 1:n
         @inbounds z = v[i] - m
         @inbounds s += abs2(z * z) * wv[i]
@@ -284,7 +284,7 @@ end
 
 function _momentk(v::RealArray, k::Int, m::Real)
     n = length(v)
-    s = 0.0
+    s = zero(eltype(v))
     for i = 1:n
         @inbounds z = v[i] - m
         s += (z ^ k)
@@ -294,7 +294,7 @@ end
 
 function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real)
     n = length(v)
-    s = 0.0
+    s = zero(promote_type(eltype(v), eltype(wv)))
     for i = 1:n
         @inbounds z = v[i] - m
         @inbounds s += (z ^ k) * wv[i]
@@ -341,8 +341,9 @@ specifying a weighting vector `wv` and a center `m`.
 """
 function skewness(v::RealArray, m::Real)
     n = length(v)
-    cm2 = 0.0   # empirical 2nd centered moment (variance)
-    cm3 = 0.0   # empirical 3rd centered moment
+    T = eltype(v)
+    cm2 = zero(T)   # empirical 2nd centered moment (variance)
+    cm3 = zero(T)   # empirical 3rd centered moment
     for i = 1:n
         @inbounds z = v[i] - m
         z2 = z * z
@@ -358,8 +359,9 @@ end
 function skewness(v::RealArray, wv::AbstractWeights, m::Real)
     n = length(v)
     length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths."))
-    cm2 = 0.0   # empirical 2nd centered moment (variance)
-    cm3 = 0.0   # empirical 3rd centered moment
+    T = eltype(v)
+    cm2 = zero(T)   # empirical 2nd centered moment (variance)
+    cm3 = zero(T)   # empirical 3rd centered moment
 
     @inbounds for i = 1:n
         x_i = v[i]
@@ -388,8 +390,9 @@ specifying a weighting vector `wv` and a center `m`.
 """
 function kurtosis(v::RealArray, m::Real)
     n = length(v)
-    cm2 = 0.0  # empirical 2nd centered moment (variance)
-    cm4 = 0.0  # empirical 4th centered moment
+    T = eltype(v)
+    cm2 = zero(T)  # empirical 2nd centered moment (variance)
+    cm4 = zero(T)  # empirical 4th centered moment
     for i = 1:n
         @inbounds z = v[i] - m
         z2 = z * z
@@ -398,14 +401,15 @@ function kurtosis(v::RealArray, m::Real)
     end
     cm4 /= n
     cm2 /= n
-    return (cm4 / (cm2 * cm2)) - 3.0
+    return (cm4 / (cm2 * cm2)) - 3
 end
 
 function kurtosis(v::RealArray, wv::AbstractWeights, m::Real)
     n = length(v)
     length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths."))
-    cm2 = 0.0  # empirical 2nd centered moment (variance)
-    cm4 = 0.0  # empirical 4th centered moment
+    T = eltype(v)
+    cm2 = zero(T)  # empirical 2nd centered moment (variance)
+    cm4 = zero(T)  # empirical 4th centered moment
 
     @inbounds for i = 1 : n
         x_i = v[i]
@@ -419,7 +423,7 @@ function kurtosis(v::RealArray, wv::AbstractWeights, m::Real)
     sw = sum(wv)
     cm4 /= sw
     cm2 /= sw
-    return (cm4 / (cm2 * cm2)) - 3.0
+    return (cm4 / (cm2 * cm2)) - 3
 end
 
 kurtosis(v::RealArray) = kurtosis(v, mean(v))

test/moments.jl

@@ -1,4 +1,5 @@
 using StatsBase
+using Statistics
 using Test
 
 @testset "StatsBase.Moments" begin
@@ -278,4 +279,19 @@ end
     @test moment(x, 5, w) ≈ sum((x2 .- 4).^5) / 5
 end
 
+@testset "Preservation of eltypes in moments" begin
+    xs = Float16[1, 2, 3, 4, 5];
+    ws = AnalyticWeights(Float16[1, 1, 1, 1, 1]);
+    @test typeof(std(xs)) === Float16
+    @test typeof(var(xs)) === Float16
+    @test typeof(mean(xs, ws)) === Float16
+    @test typeof(std(xs, ws)) === Float16
+    @test typeof(var(xs, ws)) === Float16
+    @test typeof(skewness(xs, ws)) === Float16
+    @test typeof(kurtosis(xs, ws)) === Float16
+    for i in 1:5
+        @test typeof(moment(xs, i, ws)) === Float16
+    end
+end
+
 end # @testset "StatsBase.Moments"