Make covariance and correlation work for iterators, skipmissing in particular.

Project.toml

@@ -1,5 +1,4 @@
 name = "Statistics"
-uuid = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/Statistics.jl

@@ -494,7 +494,7 @@ unscaled_covzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int) =
     (vardim == 1 ? *(transpose(x), _conj(y)) : *(x, adjoint(y)))
 
 # covzm (with centered data)
-
+covzm(itr::Any; corrected::Bool = true) = covzm(collect(itr); corrected = corrected)
 covzm(x::AbstractVector; corrected::Bool=true) = unscaled_covzm(x) / (length(x) - Int(corrected))
 function covzm(x::AbstractMatrix, vardim::Int=1; corrected::Bool=true)
     C = unscaled_covzm(x, vardim)
@@ -504,6 +504,7 @@ function covzm(x::AbstractMatrix, vardim::Int=1; corrected::Bool=true)
     A .= A .* b
     return A
 end
+covzm(x::Any, y::Any; corrected::Bool = true) = covzm(collect(x), collect(y); corrected = corrected)
 covzm(x::AbstractVector, y::AbstractVector; corrected::Bool=true) =
     unscaled_covzm(x, y) / (length(x) - Int(corrected))
 function covzm(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int=1; corrected::Bool=true)
@@ -518,22 +519,32 @@ end
 # covm (with provided mean)
 ## Use map(t -> t - xmean, x) instead of x .- xmean to allow for Vector{Vector}
 ## which can't be handled by broadcast
+covm(itr::Any, itrmean; corrected::Bool=true) = 
+    covm(collect(itr), itrmean; corrected=corrected)
 covm(x::AbstractVector, xmean; corrected::Bool=true) =
     covzm(map(t -> t - xmean, x); corrected=corrected)
 covm(x::AbstractMatrix, xmean, vardim::Int=1; corrected::Bool=true) =
     covzm(x .- xmean, vardim; corrected=corrected)
+covm(x::Any, xmean, y::Any, ymean; corrected::Bool=true) =
+    covzm(map(t -> t - xmean, x), map(t -> t - ymean, y); corrected=corrected)
 covm(x::AbstractVector, xmean, y::AbstractVector, ymean; corrected::Bool=true) =
     covzm(map(t -> t - xmean, x), map(t -> t - ymean, y); corrected=corrected)
 covm(x::AbstractVecOrMat, xmean, y::AbstractVecOrMat, ymean, vardim::Int=1; corrected::Bool=true) =
     covzm(x .- xmean, y .- ymean, vardim; corrected=corrected)
 
 # cov (API)
 """
-    cov(x::AbstractVector; corrected::Bool=true)
+    cov(itr::Any; corrected::Bool=true)
 
-Compute the variance of the vector `x`. If `corrected` is `true` (the default) then the sum
-is scaled with `n-1`, whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x)`.
+Compute the variance of the iterator `itr`. If `corrected` is `true` (the default) then the sum
+is scaled with `n-1`, whereas the sum is scaled with `n` if `corrected` is `false` where 
+``n`` is the number of elements.
 """
+function cov(itr::Any; corrected::Bool=true)
+    x = collect(itr)
+    meanx = mean(x)
+    covzm(map!(t -> t - meanx, x, x); corrected=corrected)
+end
 cov(x::AbstractVector; corrected::Bool=true) = covm(x, mean(x); corrected=corrected)
 
 """
@@ -546,14 +557,24 @@ if `corrected` is `false` where `n = size(X, dims)`.
 cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true) =
     covm(X, _vmean(X, dims), dims; corrected=corrected)
 
+
 """
-    cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true)
+    cov(x::Any, y::Any; corrected::Bool=true)
 
-Compute the covariance between the vectors `x` and `y`. If `corrected` is `true` (the
+Compute the covariance between the iterators `x` and `y`. If `corrected` is `true` (the
 default), computes ``\\frac{1}{n-1}\\sum_{i=1}^n (x_i-\\bar x) (y_i-\\bar y)^*`` where
-``*`` denotes the complex conjugate and `n = length(x) = length(y)`. If `corrected` is
+``*`` denotes the complex conjugate and ``n`` the number of elements. If `corrected` is
 `false`, computes ``\\frac{1}{n}\\sum_{i=1}^n (x_i-\\bar x) (y_i-\\bar y)^*``.
 """
+function cov(x::Any, y::Any; corrected::Bool=true)
+    cx = collect(x)
+    cy = collect(y)
+    meanx = mean(cx)
+    meany = mean(cy)
+    dx = map!(t -> t - meanx, cx, cx)
+    dy = map!(t -> t - meany, cy, cy)
+    covzm(dx, dy; corrected=corrected)
+end
 cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true) =
     covm(x, mean(x), y, mean(y); corrected=corrected)
 
@@ -630,7 +651,13 @@ function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::AbstractArray)
 end
 
 # corzm (non-exported, with centered data)
-
+function corzm(itr::Any) 
+    if Base.IteratorEltype(itr) isa Base.HasEltype && isconcrete(eltype(itr))
+        return one(real(eltype(itr)))
+    else
+        return one(real(eltype(collect(itr))))
+    end
+end
 corzm(x::AbstractVector{T}) where {T} = one(real(T))
 function corzm(x::AbstractMatrix, vardim::Int=1)
     c = unscaled_covzm(x, vardim)
@@ -644,9 +671,10 @@ corzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int=1) =
     cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, dims=vardim)), sqrt!(sum(abs2, y, dims=vardim)))
 
 # corm
-
+corm(x::Any, xmean) = corm(collect(x), xmean)
 corm(x::AbstractVector{T}, xmean) where {T} = one(real(T))
 corm(x::AbstractMatrix, xmean, vardim::Int=1) = corzm(x .- xmean, vardim)
+corm(x::Any, mx, y::Any, my) = corm(collect(x), mx, collect(y), my)
 function corm(x::AbstractVector, mx, y::AbstractVector, my)
     require_one_based_indexing(x, y)
     n = length(x)
@@ -675,11 +703,18 @@ corm(x::AbstractVecOrMat, xmean, y::AbstractVecOrMat, ymean, vardim::Int=1) =
 
 # cor
 """
-    cor(x::AbstractVector)
+    cor(itr::Any)
 
 Return the number one.
 """
-cor(x::AbstractVector) = one(real(eltype(x)))
+function cor(itr::Any)
+    if Base.IteratorEltype(itr) isa Base.HasEltype && isconcrete(eltype(itr))
+        return one(real(eltype(itr)))
+    else
+        return one(real(eltype(collect(itr))))
+    end
+end
+cor(x::AbstractVector{T}) where {T} = one(real(T))
 
 """
     cor(X::AbstractMatrix; dims::Int=1)
@@ -688,6 +723,18 @@ Compute the Pearson correlation matrix of the matrix `X` along the dimension `di
 """
 cor(X::AbstractMatrix; dims::Int=1) = corm(X, _vmean(X, dims), dims)
 
+"""
+    cor(x::AbstractVector, y::AbstractVector)
+
+Compute the Pearson correlation between the vectors `x` and `y`.
+"""
+function cor(x::Any, y::Any)
+    cx = collect(x)
+    cy = collect(y)
+
+    corm(cx, mean(cx), cy, mean(cy))
+end
+
 """
     cor(x::AbstractVector, y::AbstractVector)
 

test/runtests.jl

@@ -339,11 +339,16 @@ Y = [6.0  2.0;
             x1 = vec(X[1,:])
             y1 = vec(Y[1,:])
         end
+        @show x1
+        x1_itr = (x1i for x1i in x1)
+        y1_itr = skipmissing(y1)
 
         c = zm ? Statistics.covm(x1, 0, corrected=cr) :
                  cov(x1, corrected=cr)
+        c_itr = zm ? Statistics.covm(x1_itr, 0, corrected=cr) :
+                     cov(x1_itr, corrected=cr)
         @test isa(c, Float64)
-        @test c ≈ Cxx[1,1]
+        @test c ≈ c_itr ≈ Cxx[1,1]
         @inferred cov(x1, corrected=cr)
 
         @test cov(X) == Statistics.covm(X, mean(X, dims=1))
@@ -356,6 +361,8 @@ Y = [6.0  2.0;
         @test cov(x1, y1) == Statistics.covm(x1, mean(x1), y1, mean(y1))
         c = zm ? Statistics.covm(x1, 0, y1, 0, corrected=cr) :
                  cov(x1, y1, corrected=cr)
+        c_itr = zm ? Statistics.covm(x1_itr, 0, y1_itr, 0, corrected=cr) :
+                     cov(x1_itr, y1_itr, corrected=cr)
         @test isa(c, Float64)
         @test c ≈ Cxy[1,1]
         @inferred cov(x1, y1, corrected=cr)
@@ -426,10 +433,13 @@ end
             x1 = vec(X[1,:])
             y1 = vec(Y[1,:])
         end
+        x1_itr = (x1i for x1i in x1)
+        y1_itr = skipmissing(y1)
 
         c = zm ? Statistics.corm(x1, 0) : cor(x1)
+        c_itr = zm ? Statistics.corm(x1_itr, 0) : cor(x1_itr)
         @test isa(c, Float64)
-        @test c ≈ Cxx[1,1]
+        @test c ≈ c_itr ≈ Cxx[1,1]
         @inferred cor(x1)
 
         @test cor(X) == Statistics.corm(X, mean(X, dims=1))
@@ -440,8 +450,9 @@ end
 
         @test cor(x1, y1) == Statistics.corm(x1, mean(x1), y1, mean(y1))
         c = zm ? Statistics.corm(x1, 0, y1, 0) : cor(x1, y1)
+        c_itr = zm ? Statistics.corm(x1_itr, 0, y1_itr, 0) : cor(x1_itr, y1_itr)
         @test isa(c, Float64)
-        @test c ≈ Cxy[1,1]
+        @test c ≈ c_itr ≈ Cxy[1,1]
         @inferred cor(x1, y1)
 
         if vd == 1