Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Make view(::AbstractWeights, ...) return an AbstractWeights #723

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions docs/src/weights.md
Original file line number Diff line number Diff line change
Expand Up @@ -8,6 +8,7 @@ In statistical applications, it is not uncommon to assign weights to samples. To
!!! note
- The weight vector is a light-weight wrapper of the input vector. The input vector is NOT copied during construction.
- The weight vector maintains the sum of weights, which is computed upon construction. If the value of the sum is pre-computed, one can supply it as the second argument to the constructor and save the time of computing the sum again.
- Views of weight vectors are also `AbstractWeights`, but their sum has to be recomputed by each call to `sum` as weights in the parent vector may have been mutated.
nalimilan marked this conversation as resolved.
Show resolved Hide resolved


## Implementations
Expand Down
55 changes: 45 additions & 10 deletions src/weights.jl
Original file line number Diff line number Diff line change
Expand Up @@ -11,14 +11,19 @@ macro weights(name)
return quote
mutable struct $name{S<:Real, T<:Real, V<:AbstractVector{T}} <: AbstractWeights{S, T, V}
values::V
sum::S
"""
Pre-computed sum. Private field, to be accessed only via `sum(wv)`.
Set to `missing` by `view` as we cannot know when the parent is mutated.
"""
sum::Union{S, Missing}
end
$(esc(name))(values::AbstractVector{<:Real}) = $(esc(name))(values, sum(values))
$(esc(:weightstype))(::Type{<:$(esc(name))}) = $(esc(name))
bkamins marked this conversation as resolved.
Show resolved Hide resolved
end
end

length(wv::AbstractWeights) = length(wv.values)
sum(wv::AbstractWeights) = wv.sum
sum(wv::AbstractWeights{S}) where {S<:Real} = wv.sum::S
isempty(wv::AbstractWeights) = isempty(wv.values)
size(wv::AbstractWeights) = size(wv.values)

Expand All @@ -35,7 +40,33 @@ end
W(v, sum(v))
end

Base.getindex(wv::W, ::Colon) where {W <: AbstractWeights} = W(copy(wv.values), sum(wv))
Base.getindex(wv::AbstractWeights, ::Colon) = copy(wv)

Base.copy(wv::W) where {W <: AbstractWeights} =
weightstype(W)(copy(wv.values), sum(wv))

@propagate_inbounds function Base.view(wv::W, inds...) where
{S <: Real, W <: AbstractWeights{S}}
@boundscheck checkbounds(wv, inds...)
@inbounds v = invoke(view, Tuple{AbstractArray, Vararg{Any}}, wv, inds...)
nalimilan marked this conversation as resolved.
Show resolved Hide resolved
weightstype(W){S, eltype(wv), typeof(v)}(v, missing)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
weightstype(W){S, eltype(wv), typeof(v)}(v, missing)
return weightstype(W){S, eltype(wv), typeof(v)}(v, missing)

end

# This method is implemented for backward compatibility
bkamins marked this conversation as resolved.
Show resolved Hide resolved
@propagate_inbounds function Base.view(wv::W, inds::Union{Integer, CartesianIndex}...) where
{S <: Real, W <: AbstractWeights{S}}
@boundscheck checkbounds(wv, inds...)
@inbounds invoke(view, Tuple{AbstractArray, Vararg{Any}}, wv, inds...)
end

# Always recompute the sum for views of AbstractWeights, as we cannot know whether
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

maybe move the definitions of sum and copy into one place to make sure the reader can see both definitions side by side?

# the parent array has been mutated
Base.sum(wv::AbstractWeights{S, T, <:SubArray}) where {S<:Real, T<:Real} =
sum(wv.values)

Base.copy(wv::W) where
{S<:Real, T<:Real, W<:AbstractWeights{S, T, <:SubArray{T, <:Any, <:AbstractWeights}}} =
weightstype(W)(copy(view(parent(wv.values).values, parentindices(wv.values)...)), sum(wv))

@propagate_inbounds function Base.setindex!(wv::AbstractWeights, v::Real, i::Int)
s = v - wv[i]
Expand Down Expand Up @@ -85,7 +116,7 @@ if `corrected=true`.
@inline function varcorrection(w::Weights, corrected::Bool=false)
corrected && throw(ArgumentError("Weights type does not support bias correction: " *
"use FrequencyWeights, AnalyticWeights or ProbabilityWeights if applicable."))
1 / w.sum
1 / sum(w)
end

@weights AnalyticWeights
Expand Down Expand Up @@ -118,7 +149,7 @@ aweights(vs::RealArray) = AnalyticWeights(vec(vs))
* `corrected=false`: ``\\frac{1}{\\sum w}``
"""
@inline function varcorrection(w::AnalyticWeights, corrected::Bool=false)
s = w.sum
s = sum(w)

if corrected
sum_sn = sum(x -> (x / s) ^ 2, w)
Expand Down Expand Up @@ -156,7 +187,7 @@ fweights(vs::RealArray) = FrequencyWeights(vec(vs))
* `corrected=false`: ``\\frac{1}{\\sum w}``
"""
@inline function varcorrection(w::FrequencyWeights, corrected::Bool=false)
s = w.sum
s = sum(w)

if corrected
1 / (s - 1)
Expand Down Expand Up @@ -194,7 +225,7 @@ pweights(vs::RealArray) = ProbabilityWeights(vec(vs))
* `corrected=false`: ``\\frac{1}{\\sum w}``
"""
@inline function varcorrection(w::ProbabilityWeights, corrected::Bool=false)
s = w.sum
s = sum(w)

if corrected
n = count(!iszero, w)
Expand Down Expand Up @@ -337,8 +368,12 @@ end

for w in (AnalyticWeights, FrequencyWeights, ProbabilityWeights, Weights)
@eval begin
Base.isequal(x::$w, y::$w) = isequal(x.sum, y.sum) && isequal(x.values, y.values)
Base.:(==)(x::$w, y::$w) = (x.sum == y.sum) && (x.values == y.values)
Base.isequal(x::$w, y::$w) =
(x.values isa SubArray || y.values isa SubArray || isequal(x.sum, y.sum)) &&
isequal(x.values, y.values)
Base.:(==)(x::$w, y::$w) =
(x.values isa SubArray || y.values isa SubArray || x.sum == y.sum) &&
nalimilan marked this conversation as resolved.
Show resolved Hide resolved
x.values == y.values
end
end

Expand Down Expand Up @@ -669,7 +704,7 @@ function quantile(v::RealVector{V}, w::AbstractWeights{W}, p::RealVector) where
isempty(p) && throw(ArgumentError("empty quantile array"))
all(x -> 0 <= x <= 1, p) || throw(ArgumentError("input probability out of [0,1] range"))

w.sum == 0 && throw(ArgumentError("weight vector cannot sum to zero"))
sum(w) == 0 && throw(ArgumentError("weight vector cannot sum to zero"))
length(v) == length(w) || throw(ArgumentError("data and weight vectors must be the same size," *
"got $(length(v)) and $(length(w))"))
for x in w.values
Expand Down
16 changes: 11 additions & 5 deletions test/cov.jl
Original file line number Diff line number Diff line change
Expand Up @@ -4,9 +4,15 @@ using LinearAlgebra, Random, Test
struct EmptyCovarianceEstimator <: CovarianceEstimator end

@testset "StatsBase.Covariance" begin

function viewweights(f)
wv -> view(f([wv; 100]), axes(wv, 1))
end

weight_funcs = (weights, aweights, fweights, pweights)

@testset "$f" for f in weight_funcs
@testset "$f with $viewf" for f in weight_funcs, viewf in (identity, viewweights)
fw = viewf(f)
X = randn(3, 8)

Z1 = X .- mean(X, dims = 1)
Expand All @@ -21,8 +27,8 @@ weight_funcs = (weights, aweights, fweights, pweights)
w2[1] += 1
end

wv1 = f(w1)
wv2 = f(w2)
wv1 = fw(w1)
wv2 = fw(w2)

Z1w = X .- mean(X, wv1, dims=1)
Z2w = X .- mean(X, wv2, dims=2)
Expand Down Expand Up @@ -237,8 +243,8 @@ weight_funcs = (weights, aweights, fweights, pweights)
end

@testset "Correlation" begin
@test cor(X, f(ones(3)), 1) ≈ cor(X, dims = 1)
@test cor(X, f(ones(8)), 2) ≈ cor(X, dims = 2)
@test cor(X, fw(ones(3)), 1) ≈ cor(X, dims = 1)
@test cor(X, fw(ones(8)), 2) ≈ cor(X, dims = 2)

cov1 = cov(X, wv1, 1; corrected=false)
std1 = std(X, wv1, 1; corrected=false)
Expand Down
39 changes: 25 additions & 14 deletions test/moments.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2,15 +2,21 @@ using StatsBase
using Test

@testset "StatsBase.Moments" begin

function viewweights(f)
wv -> view(f([wv; 100]), axes(wv, 1))
end

weight_funcs = (weights, aweights, fweights, pweights)

##### weighted var & std

x = [0.57, 0.10, 0.91, 0.72, 0.46, 0.0]
w = [3.84, 2.70, 8.29, 8.91, 9.71, 0.0]

@testset "Uncorrected with $f" for f in weight_funcs
wv = f(w)
@testset "Uncorrected with $f and $viewf" for f in weight_funcs,
viewf in (identity, viewweights)
wv = viewf(f)(w)
m = mean(x, wv)

# expected uncorrected output
Expand Down Expand Up @@ -52,8 +58,9 @@ end
expected_var = [NaN, 0.0694434191182236, 0.05466601256158146, 0.06628969012045285]
expected_std = sqrt.(expected_var)

@testset "Corrected with $(weight_funcs[i])" for i in eachindex(weight_funcs)
wv = weight_funcs[i](w)
@testset "Corrected with $(weight_funcs[i])" for i in eachindex(weight_funcs),
viewf in (identity, viewweights)
wv = viewf(weight_funcs[i])(w)
m = mean(x, wv)

@testset "Variance" begin
Expand Down Expand Up @@ -107,9 +114,10 @@ x = rand(5, 6)
w1 = [0.57, 5.10, 0.91, 1.72, 0.0]
w2 = [3.84, 2.70, 8.29, 8.91, 9.71, 0.0]

@testset "Uncorrected with $f" for f in weight_funcs
wv1 = f(w1)
wv2 = f(w2)
@testset "Uncorrected with $f and $viewf" for f in weight_funcs,
viewf in (identity, viewweights)
wv1 = viewf(f)(w1)
wv2 = viewf(f)(w2)
m1 = mean(x, wv1, dims=1)
m2 = mean(x, wv2, dims=2)

Expand Down Expand Up @@ -165,9 +173,10 @@ w2 = [3.84, 2.70, 8.29, 8.91, 9.71, 0.0]
end
end

@testset "Corrected with $f" for f in weight_funcs
wv1 = f(w1)
wv2 = f(w2)
@testset "Corrected with $f and $viewf" for f in weight_funcs,
viewf in (identity, viewweights)
wv1 = viewf(f)(w1)
wv2 = viewf(f)(w2)
m1 = mean(x, wv1, dims=1)
m2 = mean(x, wv2, dims=2)

Expand Down Expand Up @@ -241,8 +250,9 @@ end
end
end

@testset "Skewness and Kurtosis with $f" for f in weight_funcs
wv = f(ones(5) * 2.0)
@testset "Skewness and Kurtosis with $f and $viewf" for f in weight_funcs,
viewf in (identity, viewweights)
wv = viewf(f)(ones(5) * 2.0)

@test skewness(1:5) ≈ 0.0
@test skewness([1, 2, 3, 4, 5]) ≈ 0.0
Expand All @@ -258,7 +268,8 @@ end
@test kurtosis([1, 2, 3, 4, 5], wv) ≈ -1.3
end

@testset "General Moments with $f" for f in weight_funcs
@testset "General Moments with $f and $viewf" for f in weight_funcs,
viewf in (identity, viewweights)
x = collect(2.0:8.0)
@test moment(x, 2) ≈ sum((x .- 5).^2) / length(x)
@test moment(x, 3) ≈ sum((x .- 5).^3) / length(x)
Expand All @@ -270,7 +281,7 @@ end
@test moment(x, 4, 4.0) ≈ sum((x .- 4).^4) / length(x)
@test moment(x, 5, 4.0) ≈ sum((x .- 4).^5) / length(x)

w = f([1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0])
w = viewf(f)([1.0, 1.0, 1.0, 1.0, 1.0, 0.0, 0.0])
x2 = collect(2.0:6.0)
@test moment(x, 2, w) ≈ sum((x2 .- 4).^2) / 5
@test moment(x, 3, w) ≈ sum((x2 .- 4).^3) / 5
Expand Down
Loading