From fe2b61119e4e9d797c915ce6b3104a37bdc97d24 Mon Sep 17 00:00:00 2001 From: Fredrik Ekre Date: Fri, 25 May 2018 08:33:11 +0200 Subject: [PATCH] move cor, cov, std, stdm, var, varm and linreg to StatsBase fix https://github.com/JuliaLang/julia/pull/25571#discussion_r162214400 (included in https://github.com/JuliaStats/StatsBase.jl/pull/379) fix #23769 (included in https://github.com/JuliaStats/StatsBase.jl/pull/379) fix #27140 --- base/deprecated.jl | 22 +- base/exports.jl | 6 - base/statistics.jl | 481 ---------------------- base/sysimg.jl | 1 - doc/src/base/math.md | 6 - doc/src/manual/interfaces.md | 10 +- stdlib/LinearAlgebra/docs/src/index.md | 1 - stdlib/LinearAlgebra/src/LinearAlgebra.jl | 1 - stdlib/LinearAlgebra/src/generic.jl | 40 -- stdlib/LinearAlgebra/test/generic.jl | 47 --- stdlib/SparseArrays/src/linalg.jl | 24 -- stdlib/SparseArrays/src/sparsematrix.jl | 61 --- stdlib/SparseArrays/test/sparse.jl | 62 +-- test/docs.jl | 1 - test/offsetarray.jl | 3 - test/reducedim.jl | 6 - test/statistics.jl | 348 +--------------- 17 files changed, 17 insertions(+), 1103 deletions(-) diff --git a/base/deprecated.jl b/base/deprecated.jl index 13342fe6e28b7..449a80ba4de32 100644 --- a/base/deprecated.jl +++ b/base/deprecated.jl @@ -293,12 +293,6 @@ deprecate(Base, :DSP, 2) using .DSP export conv, conv2, deconv, filt, filt!, xcorr -# PR #21709 -@deprecate cov(x::AbstractVector, corrected::Bool) cov(x, corrected=corrected) -@deprecate cov(x::AbstractMatrix, vardim::Int, corrected::Bool) cov(x, dims=vardim, corrected=corrected) -@deprecate cov(X::AbstractVector, Y::AbstractVector, corrected::Bool) cov(X, Y, corrected=corrected) -@deprecate cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, vardim::Int, corrected::Bool) cov(X, Y, dims=vardim, corrected=corrected) - # PR #22325 # TODO: when this replace is removed from deprecated.jl: # 1) rename the function replace_new from strings/util.jl to replace @@ -1351,13 +1345,6 @@ export readandwrite @deprecate findmin(A::AbstractArray, dims) findmin(A, dims=dims) @deprecate mean(A::AbstractArray, dims) mean(A, dims=dims) -@deprecate varm(A::AbstractArray, m::AbstractArray, dims; kwargs...) varm(A, m; kwargs..., dims=dims) -@deprecate var(A::AbstractArray, dims; kwargs...) var(A; kwargs..., dims=dims) -@deprecate std(A::AbstractArray, dims; kwargs...) std(A; kwargs..., dims=dims) -@deprecate cov(X::AbstractMatrix, dim::Int; kwargs...) cov(X; kwargs..., dims=dim) -@deprecate cov(x::AbstractVecOrMat, y::AbstractVecOrMat, dim::Int; kwargs...) cov(x, y; kwargs..., dims=dim) -@deprecate cor(X::AbstractMatrix, dim::Int) cor(X, dims=dim) -@deprecate cor(x::AbstractVecOrMat, y::AbstractVecOrMat, dim::Int) cor(x, y, dims=dim) @deprecate median(A::AbstractArray, dims; kwargs...) median(A; kwargs..., dims=dims) @deprecate mapreducedim(f, op, A::AbstractArray, dims) mapreduce(f, op, A, dims=dims) @@ -1673,6 +1660,15 @@ end @deprecate next(s::AbstractString, i::Integer) iterate(s, i) @deprecate done(s::AbstractString, i::Integer) i > ncodeunits(s) +# #27140, #27152 +@deprecate_moved cor "StatsBase" +@deprecate_moved cov "StatsBase" +@deprecate_moved std "StatsBase" +@deprecate_moved stdm "StatsBase" +@deprecate_moved var "StatsBase" +@deprecate_moved varm "StatsBase" +@deprecate_moved linreg "StatsBase" + # issue #27093 # in src/jlfrontend.scm a call to `@deprecate` is generated for per-module `eval(m, x)` @eval Core Main.Base.@deprecate(eval(e), Core.eval(Main, e)) diff --git a/base/exports.jl b/base/exports.jl index ec29cf58bc9ec..58a9626273f08 100644 --- a/base/exports.jl +++ b/base/exports.jl @@ -633,8 +633,6 @@ export set_zero_subnormals, # statistics - cor, - cov, mean!, mean, median!, @@ -642,10 +640,6 @@ export middle, quantile!, quantile, - std, - stdm, - var, - varm, # iteration done, diff --git a/base/statistics.jl b/base/statistics.jl index a5310c3fc5d2c..7737ee5831fa5 100644 --- a/base/statistics.jl +++ b/base/statistics.jl @@ -79,487 +79,6 @@ mean(A::AbstractArray; dims=:) = _mean(A, dims) _mean(A::AbstractArray{T}, region) where {T} = mean!(reducedim_init(t -> t/2, +, A, region), A) _mean(A::AbstractArray, ::Colon) = sum(A) / _length(A) -##### variances ##### - -# faster computation of real(conj(x)*y) -realXcY(x::Real, y::Real) = x*y -realXcY(x::Complex, y::Complex) = real(x)*real(y) + imag(x)*imag(y) - -var(iterable; corrected::Bool=true, mean=nothing) = _var(iterable, corrected, mean) - -function _var(iterable, corrected::Bool, mean) - y = iterate(iterable) - if y === nothing - throw(ArgumentError("variance of empty collection undefined: $(repr(iterable))")) - end - count = 1 - value, state = y - y = iterate(iterable, state) - if mean === nothing - # Use Welford algorithm as seen in (among other places) - # Knuth's TAOCP, Vol 2, page 232, 3rd edition. - M = value / 1 - S = real(zero(M)) - while y !== nothing - value, state = y - y = iterate(iterable, state) - count += 1 - new_M = M + (value - M) / count - S = S + realXcY(value - M, value - new_M) - M = new_M - end - return S / (count - Int(corrected)) - elseif isa(mean, Number) # mean provided - # Cannot use a compensated version, e.g. the one from - # "Updating Formulae and a Pairwise Algorithm for Computing Sample Variances." - # by Chan, Golub, and LeVeque, Technical Report STAN-CS-79-773, - # Department of Computer Science, Stanford University, - # because user can provide mean value that is different to mean(iterable) - sum2 = abs2(value - mean::Number) - while y !== nothing - value, state = y - y = iterate(iterable, state) - count += 1 - sum2 += abs2(value - mean) - end - return sum2 / (count - Int(corrected)) - else - throw(ArgumentError("invalid value of mean, $(mean)::$(typeof(mean))")) - end -end - -centralizedabs2fun(m) = x -> abs2.(x - m) -centralize_sumabs2(A::AbstractArray, m) = - mapreduce(centralizedabs2fun(m), +, A) -centralize_sumabs2(A::AbstractArray, m, ifirst::Int, ilast::Int) = - mapreduce_impl(centralizedabs2fun(m), +, A, ifirst, ilast) - -function centralize_sumabs2!(R::AbstractArray{S}, A::AbstractArray, means::AbstractArray) where S - # following the implementation of _mapreducedim! at base/reducedim.jl - lsiz = check_reducedims(R,A) - isempty(R) || fill!(R, zero(S)) - isempty(A) && return R - - if has_fast_linear_indexing(A) && lsiz > 16 - nslices = div(_length(A), lsiz) - ibase = first(LinearIndices(A))-1 - for i = 1:nslices - @inbounds R[i] = centralize_sumabs2(A, means[i], ibase+1, ibase+lsiz) - ibase += lsiz - end - return R - end - indsAt, indsRt = safe_tail(axes(A)), safe_tail(axes(R)) # handle d=1 manually - keep, Idefault = Broadcast.shapeindexer(indsRt) - if reducedim1(R, A) - i1 = first(indices1(R)) - @inbounds for IA in CartesianIndices(indsAt) - IR = Broadcast.newindex(IA, keep, Idefault) - r = R[i1,IR] - m = means[i1,IR] - @simd for i in axes(A, 1) - r += abs2(A[i,IA] - m) - end - R[i1,IR] = r - end - else - @inbounds for IA in CartesianIndices(indsAt) - IR = Broadcast.newindex(IA, keep, Idefault) - @simd for i in axes(A, 1) - R[i,IR] += abs2(A[i,IA] - means[i,IR]) - end - end - end - return R -end - -function varm!(R::AbstractArray{S}, A::AbstractArray, m::AbstractArray; corrected::Bool=true) where S - if isempty(A) - fill!(R, convert(S, NaN)) - else - rn = div(_length(A), _length(R)) - Int(corrected) - centralize_sumabs2!(R, A, m) - R .= R .* (1 // rn) - end - return R -end - -""" - varm(v, m; dims, corrected::Bool=true) - -Compute the sample variance of a collection `v` with known mean(s) `m`, -optionally over the given dimensions. `m` may contain means for each dimension of -`v`. If `corrected` is `true`, then the sum is scaled with `n-1`, -whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x)`. - -!!! note - Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type - to represent missing values, and the [`skipmissing`](@ref) function to omit them. -""" -varm(A::AbstractArray, m::AbstractArray; corrected::Bool=true, dims=:) = _varm(A, m, corrected, dims) - -_varm(A::AbstractArray{T}, m, corrected::Bool, region) where {T} = - varm!(reducedim_init(t -> abs2(t)/2, +, A, region), A, m; corrected=corrected) - -varm(A::AbstractArray, m; corrected::Bool=true) = _varm(A, m, corrected, :) - -function _varm(A::AbstractArray{T}, m, corrected::Bool, ::Colon) where T - n = _length(A) - n == 0 && return typeof((abs2(zero(T)) + abs2(zero(T)))/2)(NaN) - return centralize_sumabs2(A, m) / (n - Int(corrected)) -end - - -""" - var(v; dims, corrected::Bool=true, mean=nothing) - -Compute the sample variance of a vector or array `v`, optionally along the given dimensions. -The algorithm will return an estimator of the generative distribution's variance -under the assumption that each entry of `v` is an IID drawn from that generative -distribution. This computation is equivalent to calculating `sum(abs2, v - mean(v)) / -(length(v) - 1)`. If `corrected` is `true`, then the sum is scaled with `n-1`, -whereas the sum is scaled with `n` if `corrected` is `false` where `n = length(x)`. -The mean `mean` over the region may be provided. - -!!! note - Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type - to represent missing values, and the [`skipmissing`](@ref) function to omit them. -""" -var(A::AbstractArray; corrected::Bool=true, mean=nothing, dims=:) = _var(A, corrected, mean, dims) - -_var(A::AbstractArray, corrected::Bool, mean, dims) = - varm(A, coalesce(mean, Base.mean(A, dims=dims)); corrected=corrected, dims=dims) - -_var(A::AbstractArray, corrected::Bool, mean, ::Colon) = - real(varm(A, coalesce(mean, Base.mean(A)); corrected=corrected)) - -varm(iterable, m; corrected::Bool=true) = _var(iterable, corrected, m) - -## variances over ranges - -varm(v::AbstractRange, m::AbstractArray) = range_varm(v, m) -varm(v::AbstractRange, m) = range_varm(v, m) - -function range_varm(v::AbstractRange, m) - f = first(v) - m - s = step(v) - l = length(v) - vv = f^2 * l / (l - 1) + f * s * l + s^2 * l * (2 * l - 1) / 6 - if l == 0 || l == 1 - return typeof(vv)(NaN) - end - return vv -end - -function var(v::AbstractRange) - s = step(v) - l = length(v) - vv = abs2(s) * (l + 1) * l / 12 - if l == 0 || l == 1 - return typeof(vv)(NaN) - end - return vv -end - - -##### standard deviation ##### - -function sqrt!(A::AbstractArray) - for i in eachindex(A) - @inbounds A[i] = sqrt(A[i]) - end - A -end - -stdm(A::AbstractArray, m; corrected::Bool=true) = - sqrt.(varm(A, m; corrected=corrected)) - -""" - std(v; corrected::Bool=true, mean=nothing, dims) - -Compute the sample standard deviation of a vector or array `v`, optionally along the given -dimensions. The algorithm returns an estimator of the generative distribution's standard -deviation under the assumption that each entry of `v` is an IID drawn from that generative -distribution. This computation is equivalent to calculating `sqrt(sum((v - mean(v)).^2) / -(length(v) - 1))`. A pre-computed `mean` may be provided. If `corrected` is `true`, -then the sum is scaled with `n-1`, whereas the sum is scaled with `n` if `corrected` is -`false` where `n = length(x)`. - -!!! note - Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type - to represent missing values, and the [`skipmissing`](@ref) function to omit them. -""" -std(A::AbstractArray; corrected::Bool=true, mean=nothing, dims=:) = _std(A, corrected, mean, dims) - -_std(A::AbstractArray, corrected::Bool, mean, dims) = - sqrt.(var(A; corrected=corrected, mean=mean, dims=dims)) - -_std(A::AbstractArray, corrected::Bool, mean, ::Colon) = - sqrt.(var(A; corrected=corrected, mean=mean)) - -_std(A::AbstractArray{<:AbstractFloat}, corrected::Bool, mean, dims) = - sqrt!(var(A; corrected=corrected, mean=mean, dims=dims)) - -_std(A::AbstractArray{<:AbstractFloat}, corrected::Bool, mean, ::Colon) = - sqrt.(var(A; corrected=corrected, mean=mean)) - -std(iterable; corrected::Bool=true, mean=nothing) = - sqrt(var(iterable, corrected=corrected, mean=mean)) - -""" - stdm(v, m; corrected::Bool=true) - -Compute the sample standard deviation of a vector `v` -with known mean `m`. If `corrected` is `true`, -then the sum is scaled with `n-1`, whereas the sum is -scaled with `n` if `corrected` is `false` where `n = length(x)`. - -!!! note - Julia does not ignore `NaN` values in the computation. Use the [`missing`](@ref) type - to represent missing values, and the [`skipmissing`](@ref) function to omit them. -""" -stdm(iterable, m; corrected::Bool=true) = - std(iterable, corrected=corrected, mean=m) - - -###### covariance ###### - -# auxiliary functions - -_conj(x::AbstractArray{<:Real}) = x -_conj(x::AbstractArray) = conj(x) - -_getnobs(x::AbstractVector, vardim::Int) = _length(x) -_getnobs(x::AbstractMatrix, vardim::Int) = size(x, vardim) - -function _getnobs(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int) - n = _getnobs(x, vardim) - _getnobs(y, vardim) == n || throw(DimensionMismatch("dimensions of x and y mismatch")) - return n -end - -_vmean(x::AbstractVector, vardim::Int) = mean(x) -_vmean(x::AbstractMatrix, vardim::Int) = mean(x, dims=vardim) - -# core functions - -unscaled_covzm(x::AbstractVector{<:Number}) = sum(abs2, x) -unscaled_covzm(x::AbstractVector) = sum(t -> t*t', x) -unscaled_covzm(x::AbstractMatrix, vardim::Int) = (vardim == 1 ? _conj(x'x) : x * x') - -unscaled_covzm(x::AbstractVector, y::AbstractVector) = sum(conj(y[i])*x[i] for i in eachindex(y, x)) -unscaled_covzm(x::AbstractVector, y::AbstractMatrix, vardim::Int) = - (vardim == 1 ? *(transpose(x), _conj(y)) : *(transpose(x), transpose(_conj(y)))) -unscaled_covzm(x::AbstractMatrix, y::AbstractVector, vardim::Int) = - (c = vardim == 1 ? *(transpose(x), _conj(y)) : x * _conj(y); reshape(c, length(c), 1)) -unscaled_covzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int) = - (vardim == 1 ? *(transpose(x), _conj(y)) : *(x, adjoint(y))) - -# covzm (with centered data) - -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) - T = promote_type(typeof(first(C) / 1), eltype(C)) - A = convert(AbstractMatrix{T}, C) - b = 1//(size(x, vardim) - corrected) - A .= A .* b - return A -end -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) - C = unscaled_covzm(x, y, vardim) - T = promote_type(typeof(first(C) / 1), eltype(C)) - A = convert(AbstractArray{T}, C) - b = 1//(_getnobs(x, y, vardim) - corrected) - A .= A .* b - return A -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(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::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) - -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)`. -""" -cov(x::AbstractVector; corrected::Bool=true) = covm(x, Base.mean(x); corrected=corrected) - -""" - cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true) - -Compute the covariance matrix of the matrix `X` along the dimension `dims`. 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 = 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) - -Compute the covariance between the vectors `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 -`false`, computes ``\\frac{1}{n}\\sum_{i=1}^n (x_i-\\bar x) (y_i-\\bar y)^*``. -""" -cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true) = - covm(x, Base.mean(x), y, Base.mean(y); corrected=corrected) - -""" - cov(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims::Int=1, corrected::Bool=true) - -Compute the covariance between the vectors or matrices `X` and `Y` along the dimension -`dims`. 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 = size(X, dims) = size(Y, dims)`. -""" -cov(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims::Int=1, corrected::Bool=true) = - covm(X, _vmean(X, dims), Y, _vmean(Y, dims), dims; corrected=corrected) - -##### correlation ##### - -""" - clampcor(x) - -Clamp a real correlation to between -1 and 1, leaving complex correlations unchanged -""" -clampcor(x::Real) = clamp(x, -1, 1) -clampcor(x) = x - -# cov2cor! - -function cov2cor!(C::AbstractMatrix{T}, xsd::AbstractArray) where T - nx = length(xsd) - size(C) == (nx, nx) || throw(DimensionMismatch("inconsistent dimensions")) - for j = 1:nx - for i = 1:j-1 - C[i,j] = adjoint(C[j,i]) - end - C[j,j] = oneunit(T) - for i = j+1:nx - C[i,j] = clampcor(C[i,j] / (xsd[i] * xsd[j])) - end - end - return C -end -function cov2cor!(C::AbstractMatrix, xsd, ysd::AbstractArray) - nx, ny = size(C) - length(ysd) == ny || throw(DimensionMismatch("inconsistent dimensions")) - for (j, y) in enumerate(ysd) # fixme (iter): here and in all `cov2cor!` we assume that `C` is efficiently indexed by integers - for i in 1:nx - C[i,j] = clampcor(C[i, j] / (xsd * y)) - end - end - return C -end -function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd) - nx, ny = size(C) - length(xsd) == nx || throw(DimensionMismatch("inconsistent dimensions")) - for j in 1:ny - for (i, x) in enumerate(xsd) - C[i,j] = clampcor(C[i,j] / (x * ysd)) - end - end - return C -end -function cov2cor!(C::AbstractMatrix, xsd::AbstractArray, ysd::AbstractArray) - nx, ny = size(C) - (length(xsd) == nx && length(ysd) == ny) || - throw(DimensionMismatch("inconsistent dimensions")) - for (i, x) in enumerate(xsd) - for (j, y) in enumerate(ysd) - C[i,j] = clampcor(C[i,j] / (x * y)) - end - end - return C -end - -# corzm (non-exported, with centered data) - -corzm(x::AbstractVector{T}) where {T} = one(real(T)) -function corzm(x::AbstractMatrix, vardim::Int=1) - c = unscaled_covzm(x, vardim) - return cov2cor!(c, collect(sqrt(c[i,i]) for i in 1:min(size(c)...))) -end -corzm(x::AbstractVector, y::AbstractMatrix, vardim::Int=1) = - cov2cor!(unscaled_covzm(x, y, vardim), sqrt(sum(abs2, x)), sqrt!(sum(abs2, y, dims=vardim))) -corzm(x::AbstractMatrix, y::AbstractVector, vardim::Int=1) = - cov2cor!(unscaled_covzm(x, y, vardim), sqrt!(sum(abs2, x, dims=vardim)), sqrt(sum(abs2, y))) -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::AbstractVector{T}, xmean) where {T} = one(real(T)) -corm(x::AbstractMatrix, xmean, vardim::Int=1) = corzm(x .- xmean, vardim) -function corm(x::AbstractVector, mx, y::AbstractVector, my) - n = length(x) - length(y) == n || throw(DimensionMismatch("inconsistent lengths")) - n > 0 || throw(ArgumentError("correlation only defined for non-empty vectors")) - - @inbounds begin - # Initialize the accumulators - xx = zero(sqrt(abs2(x[1]))) - yy = zero(sqrt(abs2(y[1]))) - xy = zero(x[1] * y[1]') - - @simd for i in eachindex(x, y) - xi = x[i] - mx - yi = y[i] - my - xx += abs2(xi) - yy += abs2(yi) - xy += xi * yi' - end - end - return clampcor(xy / max(xx, yy) / sqrt(min(xx, yy) / max(xx, yy))) -end - -corm(x::AbstractVecOrMat, xmean, y::AbstractVecOrMat, ymean, vardim::Int=1) = - corzm(x .- xmean, y .- ymean, vardim) - -# cor -""" - cor(x::AbstractVector) - -Return the number one. -""" -cor(x::AbstractVector) = one(real(eltype(x))) - -""" - cor(X::AbstractMatrix; dims::Int=1) - -Compute the Pearson correlation matrix of the matrix `X` along the dimension `dims`. -""" -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`. -""" -cor(x::AbstractVector, y::AbstractVector) = corm(x, Base.mean(x), y, Base.mean(y)) - -""" - cor(X::AbstractVecOrMat, Y::AbstractVecOrMat; dims=1) - -Compute the Pearson correlation between the vectors or matrices `X` and `Y` along the dimension `dims`. -""" -cor(x::AbstractVecOrMat, y::AbstractVecOrMat; dims::Int=1) = - corm(x, _vmean(x, dims), y, _vmean(y, dims), dims) - ##### median & quantiles ##### """ diff --git a/base/sysimg.jl b/base/sysimg.jl index a707be20977c2..4634082dc8191 100644 --- a/base/sysimg.jl +++ b/base/sysimg.jl @@ -747,7 +747,6 @@ end # @deprecate_stdlib kron LinearAlgebra true @deprecate_stdlib ldltfact LinearAlgebra true @deprecate_stdlib ldltfact! LinearAlgebra true - @deprecate_stdlib linreg LinearAlgebra true @deprecate_stdlib logabsdet LinearAlgebra true @deprecate_stdlib logdet LinearAlgebra true @deprecate_stdlib lu LinearAlgebra true diff --git a/doc/src/base/math.md b/doc/src/base/math.md index df67d02159969..bc1312708da11 100644 --- a/doc/src/base/math.md +++ b/doc/src/base/math.md @@ -189,15 +189,9 @@ Base.FastMath.@fastmath ```@docs Base.mean Base.mean! -Base.std -Base.stdm -Base.var -Base.varm Base.middle Base.median Base.median! Base.quantile Base.quantile! -Base.cov -Base.cor ``` diff --git a/doc/src/manual/interfaces.md b/doc/src/manual/interfaces.md index 8d4810f83e3c7..9c269aa6cff0d 100644 --- a/doc/src/manual/interfaces.md +++ b/doc/src/manual/interfaces.md @@ -65,7 +65,6 @@ julia> struct Squares end julia> Base.iterate(S::Squares, state=1) = state > S.count ? nothing : (state*state, state+1) - ``` With only [`iterate`](@ref) definition, the `Squares` type is already pretty powerful. @@ -84,17 +83,18 @@ julia> for i in Squares(7) 49 ``` -We can use many of the builtin methods that work with iterables, like [`in`](@ref), [`mean`](@ref) and [`std`](@ref): +We can use many of the builtin methods that work with iterables, like [`in`](@ref), +[`sum`](@ref) and [`mean`](@ref): ```jldoctest squaretype julia> 25 in Squares(10) true +julia> sum(Squares(100)) +338350 + julia> mean(Squares(100)) 3383.5 - -julia> std(Squares(100)) -3024.355854282583 ``` There are a few more methods we can extend to give Julia more information about this iterable diff --git a/stdlib/LinearAlgebra/docs/src/index.md b/stdlib/LinearAlgebra/docs/src/index.md index a78fdfa506e2e..1a44ef9abc837 100644 --- a/stdlib/LinearAlgebra/docs/src/index.md +++ b/stdlib/LinearAlgebra/docs/src/index.md @@ -372,7 +372,6 @@ Base.inv(::AbstractMatrix) LinearAlgebra.pinv LinearAlgebra.nullspace Base.kron -LinearAlgebra.linreg LinearAlgebra.exp(::StridedMatrix{<:LinearAlgebra.BlasFloat}) LinearAlgebra.log(::StridedMatrix) LinearAlgebra.sqrt(::StridedMatrix{<:Real}) diff --git a/stdlib/LinearAlgebra/src/LinearAlgebra.jl b/stdlib/LinearAlgebra/src/LinearAlgebra.jl index 3442861f3fc8c..61cefc7f0c161 100644 --- a/stdlib/LinearAlgebra/src/LinearAlgebra.jl +++ b/stdlib/LinearAlgebra/src/LinearAlgebra.jl @@ -103,7 +103,6 @@ export ldiv!, ldlt!, ldlt, - linreg, logabsdet, logdet, lowrankdowndate, diff --git a/stdlib/LinearAlgebra/src/generic.jl b/stdlib/LinearAlgebra/src/generic.jl index 8c22d377cb716..2c26661820766 100644 --- a/stdlib/LinearAlgebra/src/generic.jl +++ b/stdlib/LinearAlgebra/src/generic.jl @@ -1121,46 +1121,6 @@ isdiag(A::AbstractMatrix) = isbanded(A, 0, 0) isdiag(x::Number) = true -""" - linreg(x, y) - -Perform simple linear regression using Ordinary Least Squares. Returns `a` and `b` such -that `a + b*x` is the closest straight line to the given points `(x, y)`, i.e., such that -the squared error between `y` and `a + b*x` is minimized. - -# Examples -```julia -using PyPlot -x = 1.0:12.0 -y = [5.5, 6.3, 7.6, 8.8, 10.9, 11.79, 13.48, 15.02, 17.77, 20.81, 22.0, 22.99] -a, b = linreg(x, y) # Linear regression -plot(x, y, "o") # Plot (x, y) points -plot(x, a + b*x) # Plot line determined by linear regression -``` - -See also: - -`\\`, [`cov`](@ref), [`std`](@ref), [`mean`](@ref). - -""" -function linreg(x::AbstractVector, y::AbstractVector) - # Least squares given - # Y = a + b*X - # where - # b = cov(X, Y)/var(X) - # a = mean(Y) - b*mean(X) - if size(x) != size(y) - throw(DimensionMismatch("x has size $(size(x)) and y has size $(size(y)), " * - "but these must be the same size")) - end - mx = mean(x) - my = mean(y) - # don't need to worry about the scaling (n vs n - 1) since they cancel in the ratio - b = Base.covm(x, mx, y, my)/Base.varm(x, mx) - a = my - b*mx - return (a, b) -end - # BLAS-like in-place y = x*α+y function (see also the version in blas.jl # for BlasFloat Arrays) function axpy!(α, x::AbstractArray, y::AbstractArray) diff --git a/stdlib/LinearAlgebra/test/generic.jl b/stdlib/LinearAlgebra/test/generic.jl index 8bf054f7530b4..f56a95e0b9211 100644 --- a/stdlib/LinearAlgebra/test/generic.jl +++ b/stdlib/LinearAlgebra/test/generic.jl @@ -89,53 +89,6 @@ n = 5 # should be odd end end -@testset "linrange" begin - # make sure unequal input arrays throw an error - x = [2; 5; 6] - y = [3; 7; 10; 10] - @test_throws DimensionMismatch linreg(x, y) - x = [2 5 6] - y = [3; 7; 10] - @test_throws MethodError linreg(x, y) - - # check (UnitRange, Array) - x = 1:12 - y = [5.5; 6.3; 7.6; 8.8; 10.9; 11.79; 13.48; 15.02; 17.77; 20.81; 22.0; 22.99] - @test [linreg(x,y)...] ≈ [2.5559090909090867, 1.6960139860139862] - @test [linreg(view(x,1:6),view(y,1:6))...] ≈ [3.8366666666666642,1.3271428571428574] - - # check (LinRange, UnitRange) - x = range(1.0, stop=12.0, length=100) - y = -100:-1 - @test [linreg(x, y)...] ≈ [-109.0, 9.0] - - # check (UnitRange, UnitRange) - x = 1:12 - y = 12:-1:1 - @test [linreg(x, y)...] ≈ [13.0, -1.0] - - # check (LinRange, LinRange) - x = range(-5, stop=10, length=100) - y = range(50, stop=200, length=100) - @test [linreg(x, y)...] ≈ [100.0, 10.0] - - # check (Array, Array) - # Anscombe's quartet (https://en.wikipedia.org/wiki/Anscombe%27s_quartet) - x123 = [10.0; 8.0; 13.0; 9.0; 11.0; 14.0; 6.0; 4.0; 12.0; 7.0; 5.0] - y1 = [8.04; 6.95; 7.58; 8.81; 8.33; 9.96; 7.24; 4.26; 10.84; 4.82; 5.68] - @test [linreg(x123,y1)...] ≈ [3.0,0.5] atol=15e-5 - - y2 = [9.14; 8.14; 8.74; 8.77; 9.26; 8.10; 6.12; 3.10; 9.13; 7.26; 4.74] - @test [linreg(x123,y2)...] ≈ [3.0,0.5] atol=10e-3 - - y3 = [7.46; 6.77; 12.74; 7.11; 7.81; 8.84; 6.08; 5.39; 8.15; 6.42; 5.73] - @test [linreg(x123,y3)...] ≈ [3.0,0.5] atol=10e-3 - - x4 = [8.0; 8.0; 8.0; 8.0; 8.0; 8.0; 8.0; 19.0; 8.0; 8.0; 8.0] - y4 = [6.58; 5.76; 7.71; 8.84; 8.47; 7.04; 5.25; 12.50; 5.56; 7.91; 6.89] - @test [linreg(x4,y4)...] ≈ [3.0,0.5] atol=10e-3 -end - @testset "diag" begin A = Matrix(1.0I, 4, 4) @test diag(A) == fill(1, 4) diff --git a/stdlib/SparseArrays/src/linalg.jl b/stdlib/SparseArrays/src/linalg.jl index 170232f7a463c..cd2041e7198f0 100644 --- a/stdlib/SparseArrays/src/linalg.jl +++ b/stdlib/SparseArrays/src/linalg.jl @@ -1010,27 +1010,3 @@ end chol(A::SparseMatrixCSC) = error("Use cholesky() instead of chol() for sparse matrices.") eigen(A::SparseMatrixCSC) = error("Use IterativeEigensolvers.eigs() instead of eigen() for sparse matrices.") - -function Base.cov(X::SparseMatrixCSC; dims::Int=1, corrected::Bool=true) - vardim = dims - a, b = size(X) - n, p = vardim == 1 ? (a, b) : (b, a) - - # The covariance can be decomposed into two terms - # 1/(n - 1) ∑ (x_i - x̄)*(x_i - x̄)' = 1/(n - 1) (∑ x_i*x_i' - n*x̄*x̄') - # which can be evaluated via a sparse matrix-matrix product - - # Compute ∑ x_i*x_i' = X'X using sparse matrix-matrix product - out = Matrix(Base.unscaled_covzm(X, vardim)) - - # Compute x̄ - x̄ᵀ = mean(X, dims=vardim) - - # Subtract n*x̄*x̄' from X'X - @inbounds for j in 1:p, i in 1:p - out[i,j] -= x̄ᵀ[i] * x̄ᵀ[j]' * n - end - - # scale with the sample size n or the corrected sample size n - 1 - return rmul!(out, inv(n - corrected)) -end diff --git a/stdlib/SparseArrays/src/sparsematrix.jl b/stdlib/SparseArrays/src/sparsematrix.jl index 6986b12cd2cf5..47fe84a7cd194 100644 --- a/stdlib/SparseArrays/src/sparsematrix.jl +++ b/stdlib/SparseArrays/src/sparsematrix.jl @@ -3489,67 +3489,6 @@ function hash(A::SparseMatrixCSC{T}, h::UInt) where T hashrun(0, length(A)-lastidx, h) # Hash zeros at end end -## Statistics - -# This is the function that does the reduction underlying var/std -function Base.centralize_sumabs2!(R::AbstractArray{S}, A::SparseMatrixCSC{Tv,Ti}, means::AbstractArray) where {S,Tv,Ti} - lsiz = Base.check_reducedims(R,A) - size(means) == size(R) || error("size of means must match size of R") - isempty(R) || fill!(R, zero(S)) - isempty(A) && return R - - colptr = A.colptr - rowval = A.rowval - nzval = A.nzval - m = size(A, 1) - n = size(A, 2) - - if size(R, 1) == size(R, 2) == 1 - # Reduction along both columns and rows - R[1, 1] = Base.centralize_sumabs2(A, means[1]) - elseif size(R, 1) == 1 - # Reduction along rows - @inbounds for col = 1:n - mu = means[col] - r = convert(S, (m-colptr[col+1]+colptr[col])*abs2(mu)) - @simd for j = colptr[col]:colptr[col+1]-1 - r += abs2(nzval[j] - mu) - end - R[1, col] = r - end - elseif size(R, 2) == 1 - # Reduction along columns - rownz = fill(convert(Ti, n), m) - @inbounds for col = 1:n - @simd for j = colptr[col]:colptr[col+1]-1 - row = rowval[j] - R[row, 1] += abs2(nzval[j] - means[row]) - rownz[row] -= 1 - end - end - for i = 1:m - R[i, 1] += rownz[i]*abs2(means[i]) - end - else - # Reduction along a dimension > 2 - @inbounds for col = 1:n - lastrow = 0 - @simd for j = colptr[col]:colptr[col+1]-1 - row = rowval[j] - for i = lastrow+1:row-1 - R[i, col] = abs2(means[i, col]) - end - R[row, col] = abs2(nzval[j] - means[row, col]) - lastrow = row - end - for i = lastrow+1:m - R[i, col] = abs2(means[i, col]) - end - end - end - return R -end - ## Uniform matrix arithmetic (+)(A::SparseMatrixCSC, J::UniformScaling) = A + sparse(J, size(A)...) diff --git a/stdlib/SparseArrays/test/sparse.jl b/stdlib/SparseArrays/test/sparse.jl index fcface7889ab4..4828d732ae177 100644 --- a/stdlib/SparseArrays/test/sparse.jl +++ b/stdlib/SparseArrays/test/sparse.jl @@ -489,7 +489,7 @@ end pA = sparse(rand(3, 7)) for arr in (se33, sA, pA) - for f in (sum, prod, minimum, maximum, var) + for f in (sum, prod, minimum, maximum) farr = Array(arr) @test f(arr) ≈ f(farr) @test f(arr, dims=1) ≈ f(farr, dims=1) @@ -518,9 +518,8 @@ end @test prod(sparse(Int[])) === 1 @test_throws ArgumentError minimum(sparse(Int[])) @test_throws ArgumentError maximum(sparse(Int[])) - @test var(sparse(Int[])) === NaN - for f in (sum, prod, var) + for f in (sum, prod) @test isequal(f(spzeros(0, 1), dims=1), f(Matrix{Int}(I, 0, 1), dims=1)) @test isequal(f(spzeros(0, 1), dims=2), f(Matrix{Int}(I, 0, 1), dims=2)) @test isequal(f(spzeros(0, 1), dims=(1, 2)), f(Matrix{Int}(I, 0, 1), dims=(1, 2))) @@ -2033,63 +2032,6 @@ end @test issymmetric(B) end -# Faster covariance function for sparse matrices -# Prevents densifying the input matrix when subtracting the mean -# Test against dense implementation -# PR https://github.com/JuliaLang/julia/pull/22735 -# Part of this test needed to be hacked due to the treatment -# of Inf in sparse matrix algebra -# https://github.com/JuliaLang/julia/issues/22921 -# The issue will be resolved in -# https://github.com/JuliaLang/julia/issues/22733 -@testset "optimizing sparse $elty covariance" for elty in (Float64, Complex{Float64}) - n = 10 - p = 5 - np2 = div(n*p, 2) - nzvals, x_sparse = guardsrand(1) do - if elty <: Real - nzvals = randn(np2) - else - nzvals = complex.(randn(np2), randn(np2)) - end - nzvals, sparse(rand(1:n, np2), rand(1:p, np2), nzvals, n, p) - end - x_dense = convert(Matrix{elty}, x_sparse) - @testset "Test with no Infs and NaNs, vardim=$vardim, corrected=$corrected" for vardim in (1, 2), - corrected in (true, false) - @test cov(x_sparse, dims=vardim, corrected=corrected) ≈ - cov(x_dense , dims=vardim, corrected=corrected) - end - - @testset "Test with $x11, vardim=$vardim, corrected=$corrected" for x11 in (NaN, Inf), - vardim in (1, 2), - corrected in (true, false) - x_sparse[1,1] = x11 - x_dense[1 ,1] = x11 - - cov_sparse = cov(x_sparse, dims=vardim, corrected=corrected) - cov_dense = cov(x_dense , dims=vardim, corrected=corrected) - @test cov_sparse[2:end, 2:end] ≈ cov_dense[2:end, 2:end] - @test isfinite.(cov_sparse) == isfinite.(cov_dense) - @test isfinite.(cov_sparse) == isfinite.(cov_dense) - end - - @testset "Test with NaN and Inf, vardim=$vardim, corrected=$corrected" for vardim in (1, 2), - corrected in (true, false) - x_sparse[1,1] = Inf - x_dense[1 ,1] = Inf - x_sparse[2,1] = NaN - x_dense[2 ,1] = NaN - - cov_sparse = cov(x_sparse, dims=vardim, corrected=corrected) - cov_dense = cov(x_dense , dims=vardim, corrected=corrected) - @test cov_sparse[(1 + vardim):end, (1 + vardim):end] ≈ - cov_dense[ (1 + vardim):end, (1 + vardim):end] - @test isfinite.(cov_sparse) == isfinite.(cov_dense) - @test isfinite.(cov_sparse) == isfinite.(cov_dense) - end -end - @testset "similar should not alias the input sparse array" begin a = sparse(rand(3,3) .+ 0.1) b = similar(a, Float32, Int32) diff --git a/test/docs.jl b/test/docs.jl index cf8bb7127f9e4..5ee2e00aacdfd 100644 --- a/test/docs.jl +++ b/test/docs.jl @@ -669,7 +669,6 @@ end @test (@repl :@r_str) !== nothing # Simple tests for apropos: -@test occursin("cor", sprint(apropos, "pearson")) @test occursin("eachindex", sprint(apropos, r"ind(exes|ices)")) using Profile @test occursin("Profile.print", sprint(apropos, "print")) diff --git a/test/offsetarray.jl b/test/offsetarray.jl index 991f366b02c75..3a334a5248e6c 100644 --- a/test/offsetarray.jl +++ b/test/offsetarray.jl @@ -374,9 +374,6 @@ I = findall(!iszero, z) @test mean(x->2x, A_3_3) == 10 @test mean(A_3_3, dims=1) == median(A_3_3, dims=1) == OffsetArray([2 5 8], (0,A_3_3.offsets[2])) @test mean(A_3_3, dims=2) == median(A_3_3, dims=2) == OffsetArray(reshape([4,5,6],(3,1)), (A_3_3.offsets[1],0)) -@test var(A_3_3) == 7.5 -@test std(A_3_3, dims=1) == OffsetArray([1 1 1], (0,A_3_3.offsets[2])) -@test std(A_3_3, dims=2) == OffsetArray(reshape([3,3,3], (3,1)), (A_3_3.offsets[1],0)) @test sum(OffsetArray(fill(1,3000), -1000)) == 3000 @test vecnorm(v) ≈ vecnorm(parent(v)) diff --git a/test/reducedim.jl b/test/reducedim.jl index 2d5e5ae4e5160..e2527e3df4fcc 100644 --- a/test/reducedim.jl +++ b/test/reducedim.jl @@ -127,7 +127,6 @@ end @test prod(A) === 1 @test_throws ArgumentError minimum(A) @test_throws ArgumentError maximum(A) - @test var(A) === NaN @test isequal(sum(A, dims=1), zeros(Int, 1, 1)) @test isequal(sum(A, dims=2), zeros(Int, 0, 1)) @@ -137,10 +136,6 @@ end @test isequal(prod(A, dims=2), fill(1, 0, 1)) @test isequal(prod(A, dims=(1, 2)), fill(1, 1, 1)) @test isequal(prod(A, dims=3), fill(1, 0, 1)) - @test isequal(var(A, dims=1), fill(NaN, 1, 1)) - @test isequal(var(A, dims=2), fill(NaN, 0, 1)) - @test isequal(var(A, dims=(1, 2)), fill(NaN, 1, 1)) - @test isequal(var(A, dims=3), fill(NaN, 0, 1)) for f in (minimum, maximum) @test_throws ArgumentError f(A, dims=1) @@ -320,7 +315,6 @@ end # issue #6672 @test sum(Real[1 2 3; 4 5.3 7.1], dims=2) == reshape([6, 16.4], 2, 1) -@test std(AbstractFloat[1,2,3], dims=1) == [1.0] @test sum(Any[1 2;3 4], dims=1) == [4 6] @test sum(Vector{Int}[[1,2],[4,3]], dims=1)[1] == [5,5] diff --git a/test/statistics.jl b/test/statistics.jl index 3d7f2c846c8f5..edc401c4c937d 100644 --- a/test/statistics.jl +++ b/test/statistics.jl @@ -80,288 +80,6 @@ end end end -@testset "var & std" begin - # edge case: empty vector - # iterable; this has to throw for type stability - @test_throws ArgumentError var(()) - @test_throws ArgumentError var((); corrected=false) - @test_throws ArgumentError var((); mean=2) - @test_throws ArgumentError var((); mean=2, corrected=false) - # reduction - @test isnan(var(Int[])) - @test isnan(var(Int[]; corrected=false)) - @test isnan(var(Int[]; mean=2)) - @test isnan(var(Int[]; mean=2, corrected=false)) - # reduction across dimensions - @test isequal(var(Int[], dims=1), [NaN]) - @test isequal(var(Int[], dims=1; corrected=false), [NaN]) - @test isequal(var(Int[], dims=1; mean=[2]), [NaN]) - @test isequal(var(Int[], dims=1; mean=[2], corrected=false), [NaN]) - - # edge case: one-element vector - # iterable - @test isnan(@inferred(var((1,)))) - @test var((1,); corrected=false) === 0.0 - @test var((1,); mean=2) === Inf - @test var((1,); mean=2, corrected=false) === 1.0 - # reduction - @test isnan(@inferred(var([1]))) - @test var([1]; corrected=false) === 0.0 - @test var([1]; mean=2) === Inf - @test var([1]; mean=2, corrected=false) === 1.0 - # reduction across dimensions - @test isequal(@inferred(var([1], dims=1)), [NaN]) - @test var([1], dims=1; corrected=false) ≈ [0.0] - @test var([1], dims=1; mean=[2]) ≈ [Inf] - @test var([1], dims=1; mean=[2], corrected=false) ≈ [1.0] - - @test var(1:8) == 6. - @test varm(1:8,1) == varm(Vector(1:8),1) - @test isnan(varm(1:1,1)) - @test isnan(var(1:1)) - @test isnan(var(1:-1)) - - @test @inferred(var(1.0:8.0)) == 6. - @test varm(1.0:8.0,1.0) == varm(Vector(1.0:8.0),1) - @test isnan(varm(1.0:1.0,1.0)) - @test isnan(var(1.0:1.0)) - @test isnan(var(1.0:-1.0)) - - @test @inferred(var(1.0f0:8.0f0)) === 6.f0 - @test varm(1.0f0:8.0f0,1.0f0) == varm(Vector(1.0f0:8.0f0),1) - @test isnan(varm(1.0f0:1.0f0,1.0f0)) - @test isnan(var(1.0f0:1.0f0)) - @test isnan(var(1.0f0:-1.0f0)) - - @test varm([1,2,3], 2) ≈ 1. - @test var([1,2,3]) ≈ 1. - @test var([1,2,3]; corrected=false) ≈ 2.0/3 - @test var([1,2,3]; mean=0) ≈ 7. - @test var([1,2,3]; mean=0, corrected=false) ≈ 14.0/3 - - @test varm((1,2,3), 2) ≈ 1. - @test var((1,2,3)) ≈ 1. - @test var((1,2,3); corrected=false) ≈ 2.0/3 - @test var((1,2,3); mean=0) ≈ 7. - @test var((1,2,3); mean=0, corrected=false) ≈ 14.0/3 - @test_throws ArgumentError var((1,2,3); mean=()) - - @test var([1 2 3 4 5; 6 7 8 9 10], dims=2) ≈ [2.5 2.5]' - @test var([1 2 3 4 5; 6 7 8 9 10], dims=2; corrected=false) ≈ [2.0 2.0]' - - @test stdm([1,2,3], 2) ≈ 1. - @test std([1,2,3]) ≈ 1. - @test std([1,2,3]; corrected=false) ≈ sqrt(2.0/3) - @test std([1,2,3]; mean=0) ≈ sqrt(7.0) - @test std([1,2,3]; mean=0, corrected=false) ≈ sqrt(14.0/3) - - @test stdm([1.0,2,3], 2) ≈ 1. - @test std([1.0,2,3]) ≈ 1. - @test std([1.0,2,3]; corrected=false) ≈ sqrt(2.0/3) - @test std([1.0,2,3]; mean=0) ≈ sqrt(7.0) - @test std([1.0,2,3]; mean=0, corrected=false) ≈ sqrt(14.0/3) - - @test std([1.0,2,3]; dims=1)[] ≈ 1. - @test std([1.0,2,3]; dims=1, corrected=false)[] ≈ sqrt(2.0/3) - @test std([1.0,2,3]; dims=1, mean=[0])[] ≈ sqrt(7.0) - @test std([1.0,2,3]; dims=1, mean=[0], corrected=false)[] ≈ sqrt(14.0/3) - - @test stdm((1,2,3), 2) ≈ 1. - @test std((1,2,3)) ≈ 1. - @test std((1,2,3); corrected=false) ≈ sqrt(2.0/3) - @test std((1,2,3); mean=0) ≈ sqrt(7.0) - @test std((1,2,3); mean=0, corrected=false) ≈ sqrt(14.0/3) - - @test std([1 2 3 4 5; 6 7 8 9 10], dims=2) ≈ sqrt.([2.5 2.5]') - @test std([1 2 3 4 5; 6 7 8 9 10], dims=2; corrected=false) ≈ sqrt.([2.0 2.0]') - - let A = ComplexF64[exp(i*im) for i in 1:10^4] - @test varm(A, 0.) ≈ sum(map(abs2, A)) / (length(A) - 1) - @test varm(A, mean(A)) ≈ var(A) - end - - @test var([1//1, 2//1]) isa Rational{Int} - @test var([1//1, 2//1], dims=1) isa Vector{Rational{Int}} - - @test std([1//1, 2//1]) isa Float64 - @test std([1//1, 2//1], dims=1) isa Vector{Float64} -end - -function safe_cov(x, y, zm::Bool, cr::Bool) - n = length(x) - if !zm - x = x .- mean(x) - y = y .- mean(y) - end - dot(vec(x), vec(y)) / (n - Int(cr)) -end -X = [1.0 5.0; - 2.0 4.0; - 3.0 6.0; - 4.0 2.0; - 5.0 1.0] -Y = [6.0 2.0; - 1.0 7.0; - 5.0 8.0; - 3.0 4.0; - 2.0 3.0] - -@testset "covariance" begin - for vd in [1, 2], zm in [true, false], cr in [true, false] - # println("vd = $vd: zm = $zm, cr = $cr") - if vd == 1 - k = size(X, 2) - Cxx = zeros(k, k) - Cxy = zeros(k, k) - for i = 1:k, j = 1:k - Cxx[i,j] = safe_cov(X[:,i], X[:,j], zm, cr) - Cxy[i,j] = safe_cov(X[:,i], Y[:,j], zm, cr) - end - x1 = vec(X[:,1]) - y1 = vec(Y[:,1]) - else - k = size(X, 1) - Cxx = zeros(k, k) - Cxy = zeros(k, k) - for i = 1:k, j = 1:k - Cxx[i,j] = safe_cov(X[i,:], X[j,:], zm, cr) - Cxy[i,j] = safe_cov(X[i,:], Y[j,:], zm, cr) - end - x1 = vec(X[1,:]) - y1 = vec(Y[1,:]) - end - - c = zm ? Base.covm(x1, 0, corrected=cr) : - cov(x1, corrected=cr) - @test isa(c, Float64) - @test c ≈ Cxx[1,1] - @inferred cov(x1, corrected=cr) - - @test cov(X) == Base.covm(X, mean(X, dims=1)) - C = zm ? Base.covm(X, 0, vd, corrected=cr) : - cov(X, dims=vd, corrected=cr) - @test size(C) == (k, k) - @test C ≈ Cxx - @inferred cov(X, dims=vd, corrected=cr) - - @test cov(x1, y1) == Base.covm(x1, mean(x1), y1, mean(y1)) - c = zm ? Base.covm(x1, 0, y1, 0, corrected=cr) : - cov(x1, y1, corrected=cr) - @test isa(c, Float64) - @test c ≈ Cxy[1,1] - @inferred cov(x1, y1, corrected=cr) - - if vd == 1 - @test cov(x1, Y) == Base.covm(x1, mean(x1), Y, mean(Y, dims=1)) - end - C = zm ? Base.covm(x1, 0, Y, 0, vd, corrected=cr) : - cov(x1, Y, dims=vd, corrected=cr) - @test size(C) == (1, k) - @test vec(C) ≈ Cxy[1,:] - @inferred cov(x1, Y, dims=vd, corrected=cr) - - if vd == 1 - @test cov(X, y1) == Base.covm(X, mean(X, dims=1), y1, mean(y1)) - end - C = zm ? Base.covm(X, 0, y1, 0, vd, corrected=cr) : - cov(X, y1, dims=vd, corrected=cr) - @test size(C) == (k, 1) - @test vec(C) ≈ Cxy[:,1] - @inferred cov(X, y1, dims=vd, corrected=cr) - - @test cov(X, Y) == Base.covm(X, mean(X, dims=1), Y, mean(Y, dims=1)) - C = zm ? Base.covm(X, 0, Y, 0, vd, corrected=cr) : - cov(X, Y, dims=vd, corrected=cr) - @test size(C) == (k, k) - @test C ≈ Cxy - @inferred cov(X, Y, dims=vd, corrected=cr) - end -end - -function safe_cor(x, y, zm::Bool) - if !zm - x = x .- mean(x) - y = y .- mean(y) - end - x = vec(x) - y = vec(y) - dot(x, y) / (sqrt(dot(x, x)) * sqrt(dot(y, y))) -end -@testset "correlation" begin - for vd in [1, 2], zm in [true, false] - # println("vd = $vd: zm = $zm") - if vd == 1 - k = size(X, 2) - Cxx = zeros(k, k) - Cxy = zeros(k, k) - for i = 1:k, j = 1:k - Cxx[i,j] = safe_cor(X[:,i], X[:,j], zm) - Cxy[i,j] = safe_cor(X[:,i], Y[:,j], zm) - end - x1 = vec(X[:,1]) - y1 = vec(Y[:,1]) - else - k = size(X, 1) - Cxx = zeros(k, k) - Cxy = zeros(k, k) - for i = 1:k, j = 1:k - Cxx[i,j] = safe_cor(X[i,:], X[j,:], zm) - Cxy[i,j] = safe_cor(X[i,:], Y[j,:], zm) - end - x1 = vec(X[1,:]) - y1 = vec(Y[1,:]) - end - - c = zm ? Base.corm(x1, 0) : cor(x1) - @test isa(c, Float64) - @test c ≈ Cxx[1,1] - @inferred cor(x1) - - @test cor(X) == Base.corm(X, mean(X, dims=1)) - C = zm ? Base.corm(X, 0, vd) : cor(X, dims=vd) - @test size(C) == (k, k) - @test C ≈ Cxx - @inferred cor(X, dims=vd) - - @test cor(x1, y1) == Base.corm(x1, mean(x1), y1, mean(y1)) - c = zm ? Base.corm(x1, 0, y1, 0) : cor(x1, y1) - @test isa(c, Float64) - @test c ≈ Cxy[1,1] - @inferred cor(x1, y1) - - if vd == 1 - @test cor(x1, Y) == Base.corm(x1, mean(x1), Y, mean(Y, dims=1)) - end - C = zm ? Base.corm(x1, 0, Y, 0, vd) : cor(x1, Y, dims=vd) - @test size(C) == (1, k) - @test vec(C) ≈ Cxy[1,:] - @inferred cor(x1, Y, dims=vd) - - if vd == 1 - @test cor(X, y1) == Base.corm(X, mean(X, dims=1), y1, mean(y1)) - end - C = zm ? Base.corm(X, 0, y1, 0, vd) : cor(X, y1, dims=vd) - @test size(C) == (k, 1) - @test vec(C) ≈ Cxy[:,1] - @inferred cor(X, y1, dims=vd) - - @test cor(X, Y) == Base.corm(X, mean(X, dims=1), Y, mean(Y, dims=1)) - C = zm ? Base.corm(X, 0, Y, 0, vd) : cor(X, Y, dims=vd) - @test size(C) == (k, k) - @test C ≈ Cxy - @inferred cor(X, Y, dims=vd) - end - - @test cor(repeat(1:17, 1, 17))[2] <= 1.0 - @test cor(1:17, 1:17) <= 1.0 - @test cor(1:17, 18:34) <= 1.0 - let tmp = range(1, stop=85, length=100) - tmp2 = Vector(tmp) - @test cor(tmp, tmp) <= 1.0 - @test cor(tmp, tmp2) <= 1.0 - end -end - @testset "quantile" begin @test quantile([1,2,3,4],0.5) == 2.5 @test quantile([1,2,3,4],[0.5]) == [2.5] @@ -383,41 +101,6 @@ let y = [0.40003674665581906, 0.4085630862624367, 0.41662034698690303, 0.4166203 @test issorted(quantile(y, range(0.01, stop=0.99, length=17))) end -@testset "variance of complex arrays (#13309)" begin - z = rand(ComplexF64, 10) - @test var(z) ≈ invoke(var, Tuple{Any}, z) ≈ cov(z) ≈ var(z,dims=1)[1] ≈ sum(abs2, z .- mean(z))/9 - @test isa(var(z), Float64) - @test isa(invoke(var, Tuple{Any}, z), Float64) - @test isa(cov(z), Float64) - @test isa(var(z,dims=1), Vector{Float64}) - @test varm(z, 0.0) ≈ invoke(varm, Tuple{Any,Float64}, z, 0.0) ≈ sum(abs2, z)/9 - @test isa(varm(z, 0.0), Float64) - @test isa(invoke(varm, Tuple{Any,Float64}, z, 0.0), Float64) - @test cor(z) === 1.0 - v = varm([1.0+2.0im], 0; corrected = false) - @test v ≈ 5 - @test isa(v, Float64) -end - -@testset "cov and cor of complex arrays (issue #21093)" begin - x = [2.7 - 3.3im, 0.9 + 5.4im, 0.1 + 0.2im, -1.7 - 5.8im, 1.1 + 1.9im] - y = [-1.7 - 1.6im, -0.2 + 6.5im, 0.8 - 10.0im, 9.1 - 3.4im, 2.7 - 5.5im] - @test cov(x, y) ≈ 4.8365 - 12.119im - @test cov(y, x) ≈ 4.8365 + 12.119im - @test cov(x, reshape(y, :, 1)) ≈ reshape([4.8365 - 12.119im], 1, 1) - @test cov(reshape(x, :, 1), y) ≈ reshape([4.8365 - 12.119im], 1, 1) - @test cov(reshape(x, :, 1), reshape(y, :, 1)) ≈ reshape([4.8365 - 12.119im], 1, 1) - @test cov([x y]) ≈ [21.779 4.8365-12.119im; - 4.8365+12.119im 54.548] - @test cor(x, y) ≈ 0.14032104449218274 - 0.35160772008699703im - @test cor(y, x) ≈ 0.14032104449218274 + 0.35160772008699703im - @test cor(x, reshape(y, :, 1)) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1) - @test cor(reshape(x, :, 1), y) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1) - @test cor(reshape(x, :, 1), reshape(y, :, 1)) ≈ reshape([0.14032104449218274 - 0.35160772008699703im], 1, 1) - @test cor([x y]) ≈ [1.0 0.14032104449218274-0.35160772008699703im - 0.14032104449218274+0.35160772008699703im 1.0] -end - @testset "Issue #17153 and PR #17154" begin a = rand(10,10) b = copy(a) @@ -429,14 +112,6 @@ end @test b == a x = mean(a, dims=2) @test b == a - x = var(a, dims=1) - @test b == a - x = var(a, dims=2) - @test b == a - x = std(a, dims=1) - @test b == a - x = std(a, dims=2) - @test b == a end # dimensional correctness @@ -448,34 +123,20 @@ using .Main.TestHelpers: Furlong @test sum(r) == sum(a) == Furlong(3) @test cumsum(r) == Furlong.([1,3]) @test mean(r) == mean(a) == median(a) == median(r) == Furlong(1.5) - @test var(r) == var(a) == Furlong{2}(0.5) - @test std(r) == std(a) == Furlong{1}(sqrt(0.5)) # Issue #21786 A = [Furlong{1}(rand(-5:5)) for i in 1:2, j in 1:2] @test mean(mean(A, dims=1), dims=2)[1] === mean(A) - @test var(A, dims=1)[1] === var(A[:, 1]) - @test std(A, dims=1)[1] === std(A[:, 1]) end # Issue #22901 -@testset "var and quantile of Any arrays" begin +@testset "quantile of Any arrays" begin x = Any[1, 2, 4, 10] y = Any[1, 2, 4, 10//1] - @test var(x) === 16.25 - @test var(y) === 65//4 - @test std(x) === sqrt(16.25) @test quantile(x, 0.5) === 3.0 @test quantile(x, 1//2) === 3//1 end -@testset "Promotion in covzm. Issue #8080" begin - A = [1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1] - @test Base.covzm(A) - mean(A, dims=1)'*mean(A, dims=1)*size(A, 1)/(size(A, 1) - 1) ≈ cov(A) - A = [1//1 -1 -1; -1 1 1; -1 1 -1; 1 -1 -1; 1 -1 1] - @test (A'A - size(A, 1)*Base.mean(A, dims=1)'*Base.mean(A, dims=1))/4 == cov(A) -end - @testset "Mean along dimension of empty array" begin a0 = zeros(0) a00 = zeros(0, 0) @@ -486,10 +147,3 @@ end @test isequal(mean(a01, dims=1) , fill(NaN, 1, 1)) @test isequal(mean(a10, dims=2) , fill(NaN, 1, 1)) end - -@testset "cov/var/std of Vector{Vector}" begin - x = [[2,4,6],[4,6,8]] - @test var(x) ≈ vec(var([x[1] x[2]], dims=2)) - @test std(x) ≈ vec(std([x[1] x[2]], dims=2)) - @test cov(x) ≈ cov([x[1] x[2]], dims=2) -end