From b54b9faa33d1c70d9f692eff4aafa64b0f1167b8 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Fri, 23 Apr 2021 19:01:50 +0200 Subject: [PATCH 1/9] Preserve eltype where possible for moments --- src/moments.jl | 74 ++++++++++++++++++++++++------------------------- test/moments.jl | 17 ++++++++++++ 2 files changed, 54 insertions(+), 37 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index c5d0ae5c0..0809be2e6 100644 --- a/src/moments.jl +++ b/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) @@ -82,17 +82,17 @@ function var!(R::AbstractArray, A::RealArray, w::AbstractWeights, dims::Int; end end -function varm(A::RealArray, w::AbstractWeights, M::RealArray, dim::Int; - corrected::DepBool=nothing) +function varm(A::RealArray{T}, w::AbstractWeights, M::RealArray, dim::Int; + corrected::DepBool=nothing) where T corrected = depcheck(:varm, corrected) - varm!(similar(A, Float64, Base.reduced_indices(axes(A), dim)), A, w, M, + varm!(similar(A, promote_type(T, 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) +function var(A::RealArray{T}, w::AbstractWeights, dim::Int; mean=nothing, + corrected::DepBool=nothing) where T corrected = depcheck(:var, corrected) - var!(similar(A, Float64, Base.reduced_indices(axes(A), dim)), A, w, dim; + var!(similar(A, promote_type(T, eltype(w)), Base.reduced_indices(axes(A), dim)), A, w, dim; mean=mean, corrected=corrected) end @@ -221,9 +221,9 @@ end ##### General central moment -function _moment2(v::RealArray, m::Real; corrected=false) +function _moment2(v::RealArray{T}, m::Real; corrected=false) where T n = length(v) - s = 0.0 + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += z * z @@ -231,9 +231,9 @@ function _moment2(v::RealArray, m::Real; corrected=false) varcorrection(n, corrected) * s end -function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) +function _moment2(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false) where T n = length(v) - s = 0.0 + s = zero(promote_type(T, eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z) * wv[i] @@ -242,9 +242,9 @@ function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) varcorrection(wv, corrected) * s end -function _moment3(v::RealArray, m::Real) +function _moment3(v::RealArray{T}, m::Real) where T n = length(v) - s = 0.0 + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += z * z * z @@ -252,9 +252,9 @@ function _moment3(v::RealArray, m::Real) s / n end -function _moment3(v::RealArray, wv::AbstractWeights, m::Real) +function _moment3(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false) where T n = length(v) - s = 0.0 + s = zero(promote_type(T, eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z * z) * wv[i] @@ -262,9 +262,9 @@ function _moment3(v::RealArray, wv::AbstractWeights, m::Real) s / sum(wv) end -function _moment4(v::RealArray, m::Real) +function _moment4(v::RealArray{T}, m::Real) where T n = length(v) - s = 0.0 + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += abs2(z * z) @@ -272,9 +272,9 @@ function _moment4(v::RealArray, m::Real) s / n end -function _moment4(v::RealArray, wv::AbstractWeights, m::Real) +function _moment4(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false) where T n = length(v) - s = 0.0 + s = zero(promote_type(T, eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += abs2(z * z) * wv[i] @@ -282,9 +282,9 @@ function _moment4(v::RealArray, wv::AbstractWeights, m::Real) s / sum(wv) end -function _momentk(v::RealArray, k::Int, m::Real) +function _momentk(v::RealArray{T}, k::Int, m::Real) where T n = length(v) - s = 0.0 + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += (z ^ k) @@ -292,9 +292,9 @@ function _momentk(v::RealArray, k::Int, m::Real) s / n end -function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real) +function _momentk(v::RealArray{T}, k::Int, wv::AbstractWeights, m::Real) where T n = length(v) - s = 0.0 + s = zero(promote_type(T, eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z ^ k) * wv[i] @@ -339,10 +339,10 @@ end Compute the standardized skewness of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function skewness(v::RealArray, m::Real) +function skewness(v::RealArray{T}, m::Real) where T n = length(v) - cm2 = 0.0 # empirical 2nd centered moment (variance) - cm3 = 0.0 # empirical 3rd centered moment + 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 @@ -355,11 +355,11 @@ function skewness(v::RealArray, m::Real) return cm3 / sqrt(cm2 * cm2 * cm2) # this is much faster than cm2^1.5 end -function skewness(v::RealArray, wv::AbstractWeights, m::Real) +function skewness(v::RealArray{T}, wv::AbstractWeights, m::Real) where T 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 + cm2 = zero(T) # empirical 2nd centered moment (variance) + cm3 = zero(T) # empirical 3rd centered moment @inbounds for i = 1:n x_i = v[i] @@ -386,10 +386,10 @@ skewness(v::RealArray, wv::AbstractWeights) = skewness(v, wv, mean(v, wv)) Compute the excess kurtosis of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function kurtosis(v::RealArray, m::Real) +function kurtosis(v::RealArray{T}, m::Real) where T n = length(v) - cm2 = 0.0 # empirical 2nd centered moment (variance) - cm4 = 0.0 # empirical 4th centered moment + 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 +398,14 @@ function kurtosis(v::RealArray, m::Real) end cm4 /= n cm2 /= n - return (cm4 / (cm2 * cm2)) - 3.0 + return (cm4 / (cm2 * cm2)) - 3 * one(T) end -function kurtosis(v::RealArray, wv::AbstractWeights, m::Real) +function kurtosis(v::RealArray{T}, wv::AbstractWeights, m::Real) where T 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 + 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 +419,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 * one(T) end kurtosis(v::RealArray) = kurtosis(v, mean(v)) diff --git a/test/moments.jl b/test/moments.jl index 97fda44ac..be24b3913 100644 --- a/test/moments.jl +++ b/test/moments.jl @@ -1,4 +1,5 @@ using StatsBase +using Statistics using Test @testset "StatsBase.Moments" begin @@ -278,4 +279,20 @@ 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" From f85a002a38469d14ce4a5ae1bc363ea288f7db3f Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Fri, 23 Apr 2021 19:03:26 +0200 Subject: [PATCH 2/9] rm extra newline --- test/moments.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/moments.jl b/test/moments.jl index be24b3913..09fc53bc2 100644 --- a/test/moments.jl +++ b/test/moments.jl @@ -294,5 +294,4 @@ end end end - end # @testset "StatsBase.Moments" From 94d26562a6b9447cdd9f4367fad63245f52d1c98 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 26 Apr 2021 12:56:44 +0200 Subject: [PATCH 3/9] simplify --- src/moments.jl | 60 +++++++++++++++++++++++++++----------------------- 1 file changed, 32 insertions(+), 28 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index 0809be2e6..0812f718f 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -82,17 +82,17 @@ function var!(R::AbstractArray, A::RealArray, w::AbstractWeights, dims::Int; end end -function varm(A::RealArray{T}, w::AbstractWeights, M::RealArray, dim::Int; - corrected::DepBool=nothing) where T +function varm(A::RealArray, w::AbstractWeights, M::RealArray, dim::Int; + corrected::DepBool=nothing) corrected = depcheck(:varm, corrected) - varm!(similar(A, promote_type(T, eltype(w)), 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{T}, w::AbstractWeights, dim::Int; mean=nothing, - corrected::DepBool=nothing) where T +function var(A::RealArray, w::AbstractWeights, dim::Int; mean=nothing, + corrected::DepBool=nothing) corrected = depcheck(:var, corrected) - var!(similar(A, promote_type(T, eltype(w)), 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 @@ -221,9 +221,9 @@ end ##### General central moment -function _moment2(v::RealArray{T}, m::Real; corrected=false) where T +function _moment2(v::RealArray, m::Real; corrected=false) n = length(v) - s = zero(T) + s = zero(eltype(v)) for i = 1:n @inbounds z = v[i] - m s += z * z @@ -231,9 +231,9 @@ function _moment2(v::RealArray{T}, m::Real; corrected=false) where T varcorrection(n, corrected) * s end -function _moment2(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false) where T +function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) n = length(v) - s = zero(promote_type(T, eltype(wv))) + s = zero(promote_type(eltype(v), eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z) * wv[i] @@ -242,9 +242,9 @@ function _moment2(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false varcorrection(wv, corrected) * s end -function _moment3(v::RealArray{T}, m::Real) where T +function _moment3(v::RealArray, m::Real) n = length(v) - s = zero(T) + s = zero(eltype(v)) for i = 1:n @inbounds z = v[i] - m s += z * z * z @@ -252,9 +252,9 @@ function _moment3(v::RealArray{T}, m::Real) where T s / n end -function _moment3(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false) where T +function _moment3(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - s = zero(promote_type(T, eltype(wv))) + s = zero(promote_type(eltype(v), eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z * z) * wv[i] @@ -262,9 +262,9 @@ function _moment3(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false s / sum(wv) end -function _moment4(v::RealArray{T}, m::Real) where T +function _moment4(v::RealArray, m::Real) n = length(v) - s = zero(T) + s = zero(eltype(v)) for i = 1:n @inbounds z = v[i] - m s += abs2(z * z) @@ -272,9 +272,9 @@ function _moment4(v::RealArray{T}, m::Real) where T s / n end -function _moment4(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false) where T +function _moment4(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) n = length(v) - s = zero(promote_type(T, eltype(wv))) + s = zero(promote_type(eltype(v), eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += abs2(z * z) * wv[i] @@ -282,9 +282,9 @@ function _moment4(v::RealArray{T}, wv::AbstractWeights, m::Real; corrected=false s / sum(wv) end -function _momentk(v::RealArray{T}, k::Int, m::Real) where T +function _momentk(v::RealArray, k::Int, m::Real) n = length(v) - s = zero(T) + s = zero(eltype(v)) for i = 1:n @inbounds z = v[i] - m s += (z ^ k) @@ -292,9 +292,9 @@ function _momentk(v::RealArray{T}, k::Int, m::Real) where T s / n end -function _momentk(v::RealArray{T}, k::Int, wv::AbstractWeights, m::Real) where T +function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real) n = length(v) - s = zero(promote_type(T, eltype(wv))) + s = zero(promote_type(eltype(v), eltype(wv))) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z ^ k) * wv[i] @@ -339,8 +339,9 @@ end Compute the standardized skewness of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function skewness(v::RealArray{T}, m::Real) where T +function skewness(v::RealArray, m::Real) n = length(v) + T = eltype(v) cm2 = zero(T) # empirical 2nd centered moment (variance) cm3 = zero(T) # empirical 3rd centered moment for i = 1:n @@ -355,9 +356,10 @@ function skewness(v::RealArray{T}, m::Real) where T return cm3 / sqrt(cm2 * cm2 * cm2) # this is much faster than cm2^1.5 end -function skewness(v::RealArray{T}, wv::AbstractWeights, m::Real) where T +function skewness(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) + T = eltype(v) cm2 = zero(T) # empirical 2nd centered moment (variance) cm3 = zero(T) # empirical 3rd centered moment @@ -386,8 +388,9 @@ skewness(v::RealArray, wv::AbstractWeights) = skewness(v, wv, mean(v, wv)) Compute the excess kurtosis of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function kurtosis(v::RealArray{T}, m::Real) where T +function kurtosis(v::RealArray, m::Real) n = length(v) + T = eltype(v) cm2 = zero(T) # empirical 2nd centered moment (variance) cm4 = zero(T) # empirical 4th centered moment for i = 1:n @@ -398,12 +401,13 @@ function kurtosis(v::RealArray{T}, m::Real) where T end cm4 /= n cm2 /= n - return (cm4 / (cm2 * cm2)) - 3 * one(T) + return (cm4 / (cm2 * cm2)) - 3 end -function kurtosis(v::RealArray{T}, wv::AbstractWeights, m::Real) where T +function kurtosis(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) + T = eltype(v) cm2 = zero(T) # empirical 2nd centered moment (variance) cm4 = zero(T) # empirical 4th centered moment @@ -419,7 +423,7 @@ function kurtosis(v::RealArray{T}, wv::AbstractWeights, m::Real) where T sw = sum(wv) cm4 /= sw cm2 /= sw - return (cm4 / (cm2 * cm2)) - 3 * one(T) + return (cm4 / (cm2 * cm2)) - 3 end kurtosis(v::RealArray) = kurtosis(v, mean(v)) From 56e0fffcbc963499e43bad20216a05e4a499367c Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 26 Apr 2021 19:46:33 +0200 Subject: [PATCH 4/9] get serious about promotion --- src/moments.jl | 40 +++++++++++++++++++++++++--------------- 1 file changed, 25 insertions(+), 15 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index 0812f718f..aff7ddf69 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -85,14 +85,16 @@ end function varm(A::RealArray, w::AbstractWeights, M::RealArray, dim::Int; corrected::DepBool=nothing) corrected = depcheck(:varm, corrected) - varm!(similar(A, promote_type(eltype(A), eltype(w)), Base.reduced_indices(axes(A), dim)), A, w, M, + T = promote_type(promote_type(eltype(A), eltype(w)), eltype(M)) + varm!(similar(A, T, 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, promote_type(eltype(A), eltype(w)), Base.reduced_indices(axes(A), dim)), A, w, dim; + T = promote_type(eltype(A), eltype(w)) + var!(similar(A, T, Base.reduced_indices(axes(A), dim)), A, w, dim; mean=mean, corrected=corrected) end @@ -223,7 +225,8 @@ end ##### General central moment function _moment2(v::RealArray, m::Real; corrected=false) n = length(v) - s = zero(eltype(v)) + T = promote_type(eltype(v), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += z * z @@ -233,7 +236,8 @@ end function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) n = length(v) - s = zero(promote_type(eltype(v), eltype(wv))) + T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z) * wv[i] @@ -244,7 +248,8 @@ end function _moment3(v::RealArray, m::Real) n = length(v) - s = zero(eltype(v)) + T = promote_type(eltype(v), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += z * z * z @@ -254,7 +259,8 @@ end function _moment3(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - s = zero(promote_type(eltype(v), eltype(wv))) + T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z * z) * wv[i] @@ -264,7 +270,8 @@ end function _moment4(v::RealArray, m::Real) n = length(v) - s = zero(eltype(v)) + T = promote_type(eltype(v), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += abs2(z * z) @@ -272,9 +279,10 @@ function _moment4(v::RealArray, m::Real) s / n end -function _moment4(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) +function _moment4(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - s = zero(promote_type(eltype(v), eltype(wv))) + T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m @inbounds s += abs2(z * z) * wv[i] @@ -284,7 +292,8 @@ end function _momentk(v::RealArray, k::Int, m::Real) n = length(v) - s = zero(eltype(v)) + T = promote_type(eltype(v), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m s += (z ^ k) @@ -294,7 +303,8 @@ end function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real) n = length(v) - s = zero(promote_type(eltype(v), eltype(wv))) + T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + s = zero(T) for i = 1:n @inbounds z = v[i] - m @inbounds s += (z ^ k) * wv[i] @@ -341,7 +351,7 @@ specifying a weighting vector `wv` and a center `m`. """ function skewness(v::RealArray, m::Real) n = length(v) - T = eltype(v) + T = promote_type(eltype(v), typeof(m)) cm2 = zero(T) # empirical 2nd centered moment (variance) cm3 = zero(T) # empirical 3rd centered moment for i = 1:n @@ -359,7 +369,7 @@ end function skewness(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - T = eltype(v) + T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) cm2 = zero(T) # empirical 2nd centered moment (variance) cm3 = zero(T) # empirical 3rd centered moment @@ -390,7 +400,7 @@ specifying a weighting vector `wv` and a center `m`. """ function kurtosis(v::RealArray, m::Real) n = length(v) - T = eltype(v) + T = promote_type(eltype(v), typeof(m)) cm2 = zero(T) # empirical 2nd centered moment (variance) cm4 = zero(T) # empirical 4th centered moment for i = 1:n @@ -407,7 +417,7 @@ end function kurtosis(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - T = eltype(v) + T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) cm2 = zero(T) # empirical 2nd centered moment (variance) cm4 = zero(T) # empirical 4th centered moment From 40e7290e1bc68dc2dc2915d994a05fa7314eee78 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 26 Apr 2021 20:07:16 +0200 Subject: [PATCH 5/9] flatten --- src/moments.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index aff7ddf69..ecb7b7b85 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -85,7 +85,7 @@ end function varm(A::RealArray, w::AbstractWeights, M::RealArray, dim::Int; corrected::DepBool=nothing) corrected = depcheck(:varm, corrected) - T = promote_type(promote_type(eltype(A), eltype(w)), eltype(M)) + T = promote_type(eltype(A), eltype(w), eltype(M)) varm!(similar(A, T, Base.reduced_indices(axes(A), dim)), A, w, M, dim; corrected=corrected) end @@ -236,7 +236,7 @@ end function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) n = length(v) - T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + T = promote_type(eltype(v), eltype(wv), typeof(m)) s = zero(T) for i = 1:n @inbounds z = v[i] - m @@ -259,7 +259,7 @@ end function _moment3(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + T = promote_type(eltype(v), eltype(wv), typeof(m)) s = zero(T) for i = 1:n @inbounds z = v[i] - m @@ -281,7 +281,7 @@ end function _moment4(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + T = promote_type(eltype(v), eltype(wv), typeof(m)) s = zero(T) for i = 1:n @inbounds z = v[i] - m @@ -303,7 +303,7 @@ end function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real) n = length(v) - T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + T = promote_type(eltype(v), eltype(wv), typeof(m)) s = zero(T) for i = 1:n @inbounds z = v[i] - m @@ -369,7 +369,7 @@ end function skewness(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + T = promote_type(eltype(v), eltype(wv), typeof(m)) cm2 = zero(T) # empirical 2nd centered moment (variance) cm3 = zero(T) # empirical 3rd centered moment @@ -417,7 +417,7 @@ end function kurtosis(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - T = promote_type(promote_type(eltype(v), eltype(wv)), typeof(m)) + T = promote_type(eltype(v), eltype(wv), typeof(m)) cm2 = zero(T) # empirical 2nd centered moment (variance) cm4 = zero(T) # empirical 4th centered moment From 338eb05578859d44594fecdb82f4f0d4109d9a91 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Mon, 26 Apr 2021 22:59:54 +0200 Subject: [PATCH 6/9] even more type stable --- src/moments.jl | 24 ++++++++---------------- 1 file changed, 8 insertions(+), 16 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index ecb7b7b85..a9500a435 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -225,8 +225,7 @@ end ##### General central moment function _moment2(v::RealArray, m::Real; corrected=false) n = length(v) - T = promote_type(eltype(v), typeof(m)) - s = zero(T) + s = zero(eltype(v)) - zero(m) for i = 1:n @inbounds z = v[i] - m s += z * z @@ -236,8 +235,7 @@ end function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) n = length(v) - T = promote_type(eltype(v), eltype(wv), typeof(m)) - s = zero(T) + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z) * wv[i] @@ -248,8 +246,7 @@ end function _moment3(v::RealArray, m::Real) n = length(v) - T = promote_type(eltype(v), typeof(m)) - s = zero(T) + s = zero(eltype(v)) - zero(m) for i = 1:n @inbounds z = v[i] - m s += z * z * z @@ -259,8 +256,7 @@ end function _moment3(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - T = promote_type(eltype(v), eltype(wv), typeof(m)) - s = zero(T) + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z * z) * wv[i] @@ -270,8 +266,7 @@ end function _moment4(v::RealArray, m::Real) n = length(v) - T = promote_type(eltype(v), typeof(m)) - s = zero(T) + s = zero(eltype(v)) - zero(m) for i = 1:n @inbounds z = v[i] - m s += abs2(z * z) @@ -281,8 +276,7 @@ end function _moment4(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - T = promote_type(eltype(v), eltype(wv), typeof(m)) - s = zero(T) + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 for i = 1:n @inbounds z = v[i] - m @inbounds s += abs2(z * z) * wv[i] @@ -292,8 +286,7 @@ end function _momentk(v::RealArray, k::Int, m::Real) n = length(v) - T = promote_type(eltype(v), typeof(m)) - s = zero(T) + s = zero(eltype(v)) - zero(m) for i = 1:n @inbounds z = v[i] - m s += (z ^ k) @@ -303,8 +296,7 @@ end function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real) n = length(v) - T = promote_type(eltype(v), eltype(wv), typeof(m)) - s = zero(T) + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 for i = 1:n @inbounds z = v[i] - m @inbounds s += (z ^ k) * wv[i] From f05e14fe43f7a138948b511b40cebe88083f7cd8 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Wed, 28 Apr 2021 08:57:09 +0000 Subject: [PATCH 7/9] Apply suggestions from code review Co-authored-by: Milan Bouchet-Valat --- src/moments.jl | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index a9500a435..6ca8bdfcf 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -225,7 +225,7 @@ end ##### General central moment function _moment2(v::RealArray, m::Real; corrected=false) n = length(v) - s = zero(eltype(v)) - zero(m) + s = (zero(eltype(v)) - zero(m))^2 for i = 1:n @inbounds z = v[i] - m s += z * z @@ -235,7 +235,7 @@ end function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) n = length(v) - s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z) * wv[i] @@ -246,7 +246,7 @@ end function _moment3(v::RealArray, m::Real) n = length(v) - s = zero(eltype(v)) - zero(m) + s = (zero(eltype(v)) - zero(m))^3 for i = 1:n @inbounds z = v[i] - m s += z * z * z @@ -256,7 +256,7 @@ end function _moment3(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^3 for i = 1:n @inbounds z = v[i] - m @inbounds s += (z * z * z) * wv[i] @@ -266,7 +266,7 @@ end function _moment4(v::RealArray, m::Real) n = length(v) - s = zero(eltype(v)) - zero(m) + s = (zero(eltype(v)) - zero(m))^4 for i = 1:n @inbounds z = v[i] - m s += abs2(z * z) @@ -276,7 +276,7 @@ end function _moment4(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) - s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^4 for i = 1:n @inbounds z = v[i] - m @inbounds s += abs2(z * z) * wv[i] @@ -296,7 +296,7 @@ end function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real) n = length(v) - s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^2 + s = zero(eltype(wv)) * (zero(eltype(v)) - zero(m))^k for i = 1:n @inbounds z = v[i] - m @inbounds s += (z ^ k) * wv[i] From eb0e7c3f92fc06dca787dc2dcf60283c5f1723ad Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Thu, 29 Apr 2021 00:08:24 +0200 Subject: [PATCH 8/9] similar type stability --- src/moments.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index 6ca8bdfcf..49c9fb72d 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -65,7 +65,7 @@ function var!(R::AbstractArray, A::RealArray, w::AbstractWeights, dims::Int; if mean == 0 varm!(R, A, w, Base.reducedim_initarray(A, dims, 0, eltype(R)), dims; corrected=corrected) - elseif mean == nothing + elseif mean === nothing varm!(R, A, w, Statistics.mean(A, w, dims=dims), dims; corrected=corrected) else # check size of mean @@ -85,7 +85,8 @@ end function varm(A::RealArray, w::AbstractWeights, M::RealArray, dim::Int; corrected::DepBool=nothing) corrected = depcheck(:varm, corrected) - T = promote_type(eltype(A), eltype(w), eltype(M)) + TA, TM, Tw = eltype(A), eltype(M), eltype(w) + T = typeof( (zero(Tw) * (zero(TA) - zero(TM)))^2 ) varm!(similar(A, T, Base.reduced_indices(axes(A), dim)), A, w, M, dim; corrected=corrected) end @@ -93,9 +94,13 @@ end function var(A::RealArray, w::AbstractWeights, dim::Int; mean=nothing, corrected::DepBool=nothing) corrected = depcheck(:var, corrected) - T = promote_type(eltype(A), eltype(w)) + # doing this here instead of in var! + # allows better type stability for the returned array + M = Statistics.mean(A, w, dims=dim) + TA, TM, Tw = eltype(A), eltype(M), eltype(w) + T = typeof( (zero(Tw) * (zero(TA) - zero(TM)))^2 ) var!(similar(A, T, Base.reduced_indices(axes(A), dim)), A, w, dim; - mean=mean, corrected=corrected) + mean=M, corrected=corrected) end ## std From aacef2a3e8d8fb4e7a3bbfd01d9da6fd0fa940f2 Mon Sep 17 00:00:00 2001 From: Phillip Alday Date: Tue, 4 May 2021 23:03:41 +0200 Subject: [PATCH 9/9] skewness and kurtosis --- src/moments.jl | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/moments.jl b/src/moments.jl index 49c9fb72d..c973d8b8d 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -348,9 +348,9 @@ specifying a weighting vector `wv` and a center `m`. """ function skewness(v::RealArray, m::Real) n = length(v) - T = promote_type(eltype(v), typeof(m)) - cm2 = zero(T) # empirical 2nd centered moment (variance) - cm3 = zero(T) # empirical 3rd centered moment + T = typeof( zero(eltype(v)) - zero(m) ) + cm2 = zero(T)^2 # empirical 2nd centered moment (variance) + cm3 = zero(T)^3 # empirical 3rd centered moment for i = 1:n @inbounds z = v[i] - m z2 = z * z @@ -366,9 +366,9 @@ end function skewness(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - T = promote_type(eltype(v), eltype(wv), typeof(m)) - cm2 = zero(T) # empirical 2nd centered moment (variance) - cm3 = zero(T) # empirical 3rd centered moment + T = typeof(zero(eltype(v)) - zero(m)) + cm2 = zero(eltype(wv)) * zero(T)^2 # empirical 2nd centered moment (variance) + cm3 = zero(eltype(wv)) * zero(T)^3 # empirical 3rd centered moment @inbounds for i = 1:n x_i = v[i] @@ -397,9 +397,9 @@ specifying a weighting vector `wv` and a center `m`. """ function kurtosis(v::RealArray, m::Real) n = length(v) - T = promote_type(eltype(v), typeof(m)) - cm2 = zero(T) # empirical 2nd centered moment (variance) - cm4 = zero(T) # empirical 4th centered moment + T = typeof( zero(eltype(v)) - zero(m) ) + cm2 = zero(T)^2 # empirical 2nd centered moment (variance) + cm4 = zero(T)^4 # empirical 4th centered moment for i = 1:n @inbounds z = v[i] - m z2 = z * z @@ -414,9 +414,9 @@ end function kurtosis(v::RealArray, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) - T = promote_type(eltype(v), eltype(wv), typeof(m)) - cm2 = zero(T) # empirical 2nd centered moment (variance) - cm4 = zero(T) # empirical 4th centered moment + T = typeof(zero(eltype(v)) - zero(m)) + cm2 = zero(eltype(wv)) * zero(T)^2 # empirical 2nd centered moment (variance) + cm4 = zero(eltype(wv)) * zero(T)^4 # empirical 4th centered moment @inbounds for i = 1 : n x_i = v[i]