From 32a14aec88bfefce9495e29a2c9aa8ed6426c6f6 Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Tue, 6 Dec 2022 17:48:03 +0100 Subject: [PATCH 01/15] Get rid of old type aliases (#840) --- docs/src/counts.md | 2 +- docs/src/scalarstats.md | 2 +- src/common.jl | 25 +------------ src/counts.jl | 74 +++++++++++++++++++------------------- src/cov.jl | 6 ++-- src/deprecates.jl | 22 ++++++------ src/empirical.jl | 4 +-- src/misc.jl | 9 +++-- src/moments.jl | 78 ++++++++++++++++++++--------------------- src/rankcorr.jl | 24 ++++++------- src/ranking.jl | 8 ++--- src/sampling.jl | 20 +++++------ src/signalcorr.jl | 70 ++++++++++++++++++------------------ src/weights.jl | 34 +++++++++--------- test/rankcorr.jl | 14 ++++---- 15 files changed, 183 insertions(+), 209 deletions(-) diff --git a/docs/src/counts.md b/docs/src/counts.md index 604f7926a..e4498a1d7 100644 --- a/docs/src/counts.md +++ b/docs/src/counts.md @@ -7,7 +7,7 @@ The package provides functions to count the occurrences of distinct values. ```@docs counts proportions -addcounts!(r::AbstractArray, x::StatsBase.IntegerArray, levels::StatsBase.IntUnitRange) +addcounts!(r::AbstractArray, x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) ``` ## Counting over arbitrary distinct values diff --git a/docs/src/scalarstats.md b/docs/src/scalarstats.md index 98a73c1b1..406e4765a 100644 --- a/docs/src/scalarstats.md +++ b/docs/src/scalarstats.md @@ -69,7 +69,7 @@ percentile iqr nquantile quantile -Statistics.median(v::StatsBase.RealVector, w::AbstractWeights{<:Real}) +Statistics.median(v::AbstractVector{<:Real}, w::AbstractWeights{<:Real}) quantilerank percentilerank ``` diff --git a/src/common.jl b/src/common.jl index 17ee3a890..e4d4f712b 100644 --- a/src/common.jl +++ b/src/common.jl @@ -1,29 +1,6 @@ # common utilities -## convenient type alias -# -# These types signficantly reduces the need of using -# type parameters in functions (which are often just -# for the purpose of restricting the arrays to real) -# -# These could be removed when the Base supports -# covariant type notation, i.e. AbstractVector{<:Real} -# - -const RealArray{T<:Real,N} = AbstractArray{T,N} -const RealVector{T<:Real} = AbstractArray{T,1} -const RealMatrix{T<:Real} = AbstractArray{T,2} - -const IntegerArray{T<:Integer,N} = AbstractArray{T,N} -const IntegerVector{T<:Integer} = AbstractArray{T,1} -const IntegerMatrix{T<:Integer} = AbstractArray{T,2} - -const RealFP = Union{Float32, Float64} - -# A convenient typealias for deprecating default corrected Bool -const DepBool = Union{Bool, Nothing} - -function depcheck(fname::Symbol, varname::Symbol, b::DepBool) +function depcheck(fname::Symbol, varname::Symbol, b::Union{Bool, Nothing}) if b === nothing msg = "$fname will default to $varname=true in the future. Use $varname=false for previous behaviour." Base.depwarn(msg, fname) diff --git a/src/counts.jl b/src/counts.jl index 489196f68..d7485fd30 100644 --- a/src/counts.jl +++ b/src/counts.jl @@ -6,8 +6,6 @@ # ################################################# -const IntUnitRange{T<:Integer} = UnitRange{T} - if isdefined(Base, :ht_keyindex2) const ht_keyindex2! = Base.ht_keyindex2 else @@ -24,14 +22,14 @@ array `r`. For each `xi ∈ x`, if `xi == levels[j]`, then we increment `r[j]`. If a weighting vector `wv` is specified, the sum of weights is used rather than the raw counts. """ -function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntUnitRange) +function addcounts!(r::AbstractArray, x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) # add counts of integers from x that fall within levels to r checkbounds(r, axes(levels)...) m0 = first(levels) m1 = last(levels) - b = m0 - firstindex(levels) # firstindex(levels) == 1 because levels::IntUnitRange + b = m0 - firstindex(levels) # firstindex(levels) == 1 because levels::UnitRange{<:Integer} @inbounds for xi in x if m0 <= xi <= m1 @@ -41,7 +39,7 @@ function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntUnitRange) return r end -function addcounts!(r::AbstractArray, x::IntegerArray, levels::IntUnitRange, wv::AbstractWeights) +function addcounts!(r::AbstractArray, x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}, wv::AbstractWeights) # add wv weighted counts of integers from x that fall within levels to r length(x) == length(wv) || @@ -82,14 +80,14 @@ The output is a vector of length `length(levels)`. """ function counts end -counts(x::IntegerArray, levels::IntUnitRange) = +counts(x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) = addcounts!(zeros(Int, length(levels)), x, levels) -counts(x::IntegerArray, levels::IntUnitRange, wv::AbstractWeights) = +counts(x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}, wv::AbstractWeights) = addcounts!(zeros(eltype(wv), length(levels)), x, levels, wv) -counts(x::IntegerArray, k::Integer) = counts(x, 1:k) -counts(x::IntegerArray, k::Integer, wv::AbstractWeights) = counts(x, 1:k, wv) -counts(x::IntegerArray) = counts(x, span(x)) -counts(x::IntegerArray, wv::AbstractWeights) = counts(x, span(x), wv) +counts(x::AbstractArray{<:Integer}, k::Integer) = counts(x, 1:k) +counts(x::AbstractArray{<:Integer}, k::Integer, wv::AbstractWeights) = counts(x, 1:k, wv) +counts(x::AbstractArray{<:Integer}) = counts(x, span(x)) +counts(x::AbstractArray{<:Integer}, wv::AbstractWeights) = counts(x, span(x), wv) """ @@ -101,8 +99,8 @@ Equivalent to `counts(x, levels) / length(x)`. If a vector of weights `wv` is provided, the proportion of weights is computed rather than the proportion of raw counts. """ -proportions(x::IntegerArray, levels::IntUnitRange) = counts(x, levels) .* inv(length(x)) -proportions(x::IntegerArray, levels::IntUnitRange, wv::AbstractWeights) = +proportions(x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) = counts(x, levels) .* inv(length(x)) +proportions(x::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}, wv::AbstractWeights) = counts(x, levels, wv) .* inv(sum(wv)) """ @@ -113,14 +111,14 @@ Return the proportion of integers in 1 to `k` that occur in `x`. If a vector of weights `wv` is provided, the proportion of weights is computed rather than the proportion of raw counts. """ -proportions(x::IntegerArray, k::Integer) = proportions(x, 1:k) -proportions(x::IntegerArray, k::Integer, wv::AbstractWeights) = proportions(x, 1:k, wv) -proportions(x::IntegerArray) = proportions(x, span(x)) -proportions(x::IntegerArray, wv::AbstractWeights) = proportions(x, span(x), wv) +proportions(x::AbstractArray{<:Integer}, k::Integer) = proportions(x, 1:k) +proportions(x::AbstractArray{<:Integer}, k::Integer, wv::AbstractWeights) = proportions(x, 1:k, wv) +proportions(x::AbstractArray{<:Integer}) = proportions(x, span(x)) +proportions(x::AbstractArray{<:Integer}, wv::AbstractWeights) = proportions(x, span(x), wv) #### functions for counting a single list of integers (2D) -function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}) +function addcounts!(r::AbstractArray, x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, levels::NTuple{2,UnitRange{<:Integer}}) # add counts of pairs from zip(x,y) to r xlevels, ylevels = levels @@ -146,8 +144,8 @@ function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray, levels:: return r end -function addcounts!(r::AbstractArray, x::IntegerArray, y::IntegerArray, - levels::NTuple{2,IntUnitRange}, wv::AbstractWeights) +function addcounts!(r::AbstractArray, x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, + levels::NTuple{2,UnitRange{<:Integer}}, wv::AbstractWeights) # add counts of pairs from zip(x,y) to r length(x) == length(y) == length(wv) || @@ -182,43 +180,43 @@ end # facet functions -function counts(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}) +function counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, levels::NTuple{2,UnitRange{<:Integer}}) addcounts!(zeros(Int, length(levels[1]), length(levels[2])), x, y, levels) end -function counts(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}, wv::AbstractWeights) +function counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, levels::NTuple{2,UnitRange{<:Integer}}, wv::AbstractWeights) addcounts!(zeros(eltype(wv), length(levels[1]), length(levels[2])), x, y, levels, wv) end -counts(x::IntegerArray, y::IntegerArray, levels::IntUnitRange) = +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}) = counts(x, y, (levels, levels)) -counts(x::IntegerArray, y::IntegerArray, levels::IntUnitRange, wv::AbstractWeights) = +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, levels::UnitRange{<:Integer}, wv::AbstractWeights) = counts(x, y, (levels, levels), wv) -counts(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}) = +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, ks::NTuple{2,Integer}) = counts(x, y, (1:ks[1], 1:ks[2])) -counts(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}, wv::AbstractWeights) = +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, ks::NTuple{2,Integer}, wv::AbstractWeights) = counts(x, y, (1:ks[1], 1:ks[2]), wv) -counts(x::IntegerArray, y::IntegerArray, k::Integer) = counts(x, y, (1:k, 1:k)) -counts(x::IntegerArray, y::IntegerArray, k::Integer, wv::AbstractWeights) = +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, k::Integer) = counts(x, y, (1:k, 1:k)) +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, k::Integer, wv::AbstractWeights) = counts(x, y, (1:k, 1:k), wv) -counts(x::IntegerArray, y::IntegerArray) = counts(x, y, (span(x), span(y))) -counts(x::IntegerArray, y::IntegerArray, wv::AbstractWeights) = counts(x, y, (span(x), span(y)), wv) +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}) = counts(x, y, (span(x), span(y))) +counts(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, wv::AbstractWeights) = counts(x, y, (span(x), span(y)), wv) -proportions(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}) = +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, levels::NTuple{2,UnitRange{<:Integer}}) = counts(x, y, levels) .* inv(length(x)) -proportions(x::IntegerArray, y::IntegerArray, levels::NTuple{2,IntUnitRange}, wv::AbstractWeights) = +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, levels::NTuple{2,UnitRange{<:Integer}}, wv::AbstractWeights) = counts(x, y, levels, wv) .* inv(sum(wv)) -proportions(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}) = +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, ks::NTuple{2,Integer}) = proportions(x, y, (1:ks[1], 1:ks[2])) -proportions(x::IntegerArray, y::IntegerArray, ks::NTuple{2,Integer}, wv::AbstractWeights) = +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, ks::NTuple{2,Integer}, wv::AbstractWeights) = proportions(x, y, (1:ks[1], 1:ks[2]), wv) -proportions(x::IntegerArray, y::IntegerArray, k::Integer) = proportions(x, y, (1:k, 1:k)) -proportions(x::IntegerArray, y::IntegerArray, k::Integer, wv::AbstractWeights) = +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, k::Integer) = proportions(x, y, (1:k, 1:k)) +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, k::Integer, wv::AbstractWeights) = proportions(x, y, (1:k, 1:k), wv) -proportions(x::IntegerArray, y::IntegerArray) = proportions(x, y, (span(x), span(y))) -proportions(x::IntegerArray, y::IntegerArray, wv::AbstractWeights) = +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}) = proportions(x, y, (span(x), span(y))) +proportions(x::AbstractArray{<:Integer}, y::AbstractArray{<:Integer}, wv::AbstractWeights) = proportions(x, y, (span(x), span(y)), wv) diff --git a/src/cov.jl b/src/cov.jl index 040d3098e..a8b3cbcb8 100644 --- a/src/cov.jl +++ b/src/cov.jl @@ -91,11 +91,11 @@ _scattermatm(x::DenseMatrix, wv::AbstractWeights, mean, dims::Int) = ## weighted cov covm(x::DenseMatrix, mean, w::AbstractWeights, dims::Int=1; - corrected::DepBool=nothing) = + corrected::Union{Bool, Nothing}=nothing) = rmul!(scattermat(x, w, mean=mean, dims=dims), varcorrection(w, depcheck(:covm, :corrected, corrected))) -cov(x::DenseMatrix, w::AbstractWeights, dims::Int=1; corrected::DepBool=nothing) = +cov(x::DenseMatrix, w::AbstractWeights, dims::Int=1; corrected::Union{Bool, Nothing}=nothing) = covm(x, mean(x, w, dims=dims), w, dims; corrected=depcheck(:cov, :corrected, corrected)) function corm(x::DenseMatrix, mean, w::AbstractWeights, vardim::Int=1) @@ -118,7 +118,7 @@ function mean_and_cov(x::DenseMatrix, dims::Int=1; corrected::Bool=true) return m, covm(x, m, dims, corrected=corrected) end function mean_and_cov(x::DenseMatrix, wv::AbstractWeights, dims::Int=1; - corrected::DepBool=nothing) + corrected::Union{Bool, Nothing}=nothing) m = mean(x, wv, dims=dims) return m, cov(x, wv, dims; corrected=depcheck(:mean_and_cov, :corrected, corrected)) end diff --git a/src/deprecates.jl b/src/deprecates.jl index 9d0cc0e6e..309b93dd0 100644 --- a/src/deprecates.jl +++ b/src/deprecates.jl @@ -27,12 +27,12 @@ end @deprecate mean!(R::AbstractArray, A::AbstractArray, w::AbstractWeights, dims::Int) mean!(R, A, w, dims=dims) @deprecate mean(A::AbstractArray{T}, w::AbstractWeights{W}, dims::Int) where {T<:Number,W<:Real} mean(A, w, dims=dims) -@deprecate wquantile(v::RealVector, w::AbstractWeights{<:Real}, p::RealVector) quantile(v, w, p) -@deprecate wquantile(v::RealVector, w::AbstractWeights{<:Real}, p::Number) quantile(v, w, [p])[1] -@deprecate wquantile(v::RealVector, w::RealVector, p::RealVector) quantile(v, pweights(w), p) -@deprecate wquantile(v::RealVector, w::RealVector, p::Number) quantile(v, pweights(w), [p])[1] -@deprecate wmedian(v::RealVector, w::AbstractWeights{<:Real}) median(v, w) -@deprecate wmedian(v::RealVector, w::RealVector) median(v, weights(w)) +@deprecate wquantile(v::AbstractVector{<:Real}, w::AbstractWeights{<:Real}, p::AbstractVector{<:Real}) quantile(v, w, p) +@deprecate wquantile(v::AbstractVector{<:Real}, w::AbstractWeights{<:Real}, p::Number) quantile(v, w, [p])[1] +@deprecate wquantile(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}, p::AbstractVector{<:Real}) quantile(v, pweights(w), p) +@deprecate wquantile(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}, p::Number) quantile(v, pweights(w), [p])[1] +@deprecate wmedian(v::AbstractVector{<:Real}, w::AbstractWeights{<:Real}) median(v, w) +@deprecate wmedian(v::AbstractVector{<:Real}, w::AbstractVector{<:Real}) median(v, weights(w)) @deprecate quantile(v::AbstractArray{<:Real}) quantile(v, [.0, .25, .5, .75, 1.0]) @@ -41,8 +41,8 @@ end @deprecate values(wv::AbstractWeights) convert(Vector, wv) ### Deprecated November 2021 -@deprecate stdm(x::RealArray, w::AbstractWeights, m::Real; corrected::DepBool=nothing) std(x, w, mean=m, corrected=corrected) false -@deprecate varm(x::RealArray, w::AbstractWeights, m::Real; corrected::DepBool=nothing) var(x, w, mean=m, corrected=corrected) false -@deprecate stdm(x::RealArray, w::AbstractWeights, m::RealArray, dim::Int; corrected::DepBool=nothing) std(x, w, dim, mean=m, corrected=corrected) false -@deprecate varm(x::RealArray, w::AbstractWeights, m::RealArray, dim::Int; corrected::DepBool=nothing) var(x, w, dim, mean=m, corrected=corrected) false -@deprecate varm!(R::AbstractArray, x::RealArray, w::AbstractWeights, m::RealArray, dim::Int; corrected::DepBool=nothing) var!(R, x, w, dim, mean=m, corrected=corrected) false \ No newline at end of file +@deprecate stdm(x::AbstractArray{<:Real}, w::AbstractWeights, m::Real; corrected::Union{Bool, Nothing}=nothing) std(x, w, mean=m, corrected=corrected) false +@deprecate varm(x::AbstractArray{<:Real}, w::AbstractWeights, m::Real; corrected::Union{Bool, Nothing}=nothing) var(x, w, mean=m, corrected=corrected) false +@deprecate stdm(x::AbstractArray{<:Real}, w::AbstractWeights, m::AbstractArray{<:Real}, dim::Int; corrected::Union{Bool, Nothing}=nothing) std(x, w, dim, mean=m, corrected=corrected) false +@deprecate varm(x::AbstractArray{<:Real}, w::AbstractWeights, m::AbstractArray{<:Real}, dim::Int; corrected::Union{Bool, Nothing}=nothing) var(x, w, dim, mean=m, corrected=corrected) false +@deprecate varm!(R::AbstractArray, x::AbstractArray{<:Real}, w::AbstractWeights, m::AbstractArray{<:Real}, dim::Int; corrected::Union{Bool, Nothing}=nothing) var!(R, x, w, dim, mean=m, corrected=corrected) false diff --git a/src/empirical.jl b/src/empirical.jl index 98ef7d91b..45f985468 100644 --- a/src/empirical.jl +++ b/src/empirical.jl @@ -16,7 +16,7 @@ function (ecdf::ECDF)(x::Real) partialsum / weightsum end -function (ecdf::ECDF)(v::RealVector) +function (ecdf::ECDF)(v::AbstractVector{<:Real}) evenweights = isempty(ecdf.weights) weightsum = evenweights ? length(ecdf.sorted_values) : sum(ecdf.weights) ord = sortperm(v) @@ -53,7 +53,7 @@ evaluate CDF values on other samples. `extrema`, `minimum`, and `maximum` are supported to for obtaining the range over which function is inside the interval ``(0,1)``; the function is defined for the whole real line. """ -function ecdf(X::RealVector; weights::AbstractVector{<:Real}=Weights(Float64[])) +function ecdf(X::AbstractVector{<:Real}; weights::AbstractVector{<:Real}=Weights(Float64[])) any(isnan, X) && throw(ArgumentError("ecdf can not include NaN values")) isempty(weights) || length(X) == length(weights) || throw(ArgumentError("data and weight vectors must be the same size," * "got $(length(X)) and $(length(weights))")) diff --git a/src/misc.jl b/src/misc.jl index 25fc0f3d4..cd661c4ee 100644 --- a/src/misc.jl +++ b/src/misc.jl @@ -57,7 +57,7 @@ Reconstruct a vector from its run-length encoding (see [`rle`](@ref)). `vals` is a vector of the values and `lens` is a vector of the corresponding run lengths. """ -function inverse_rle(vals::AbstractVector{T}, lens::IntegerVector) where T +function inverse_rle(vals::AbstractVector{T}, lens::AbstractVector{<:Integer}) where T m = length(vals) length(lens) == m || raise_dimerror() @@ -131,7 +131,7 @@ julia> indicatormat([1 2 2], 2) 0 1 1 ``` """ -function indicatormat(x::IntegerArray, k::Integer; sparse::Bool=false) +function indicatormat(x::AbstractArray{<:Integer}, k::Integer; sparse::Bool=false) sparse ? _indicatormat_sparse(x, k) : _indicatormat_dense(x, k) end @@ -151,7 +151,7 @@ indicatormat(x::AbstractArray; sparse::Bool=false) = indicatormat(x, sort!(unique(x)); sparse=sparse) -function _indicatormat_dense(x::IntegerArray, k::Integer) +function _indicatormat_dense(x::AbstractArray{<:Integer}, k::Integer) n = length(x) r = zeros(Bool, k, n) for i = 1 : n @@ -174,7 +174,7 @@ function _indicatormat_dense(x::AbstractArray{T}, c::AbstractArray{T}) where T return r end -_indicatormat_sparse(x::IntegerArray, k::Integer) = (n = length(x); sparse(x, 1:n, true, k, n)) +_indicatormat_sparse(x::AbstractArray{<:Integer}, k::Integer) = (n = length(x); sparse(x, 1:n, true, k, n)) function _indicatormat_sparse(x::AbstractArray{T}, c::AbstractArray{T}) where T d = indexmap(c) @@ -187,4 +187,3 @@ function _indicatormat_sparse(x::AbstractArray{T}, c::AbstractArray{T}) where T end return sparse(rinds, 1:n, true, m, n) end - diff --git a/src/moments.jl b/src/moments.jl index b67acfb04..02988a94c 100644 --- a/src/moments.jl +++ b/src/moments.jl @@ -18,8 +18,8 @@ replacing ``\\frac{1}{\\sum{w}}`` with a factor dependent on the type of weights * `ProbabilityWeights`: ``\\frac{n}{(n - 1) \\sum w}`` where ``n`` equals `count(!iszero, w)` * `Weights`: `ArgumentError` (bias correction not supported) """ -function var(v::RealArray, w::AbstractWeights; mean=nothing, - corrected::DepBool=nothing) +function var(v::AbstractArray{<:Real}, w::AbstractWeights; mean=nothing, + corrected::Union{Bool, Nothing}=nothing) corrected = depcheck(:var, :corrected, corrected) if mean == nothing @@ -31,8 +31,8 @@ end ## var along dim -function var!(R::AbstractArray, A::RealArray, w::AbstractWeights, dims::Int; - mean=nothing, corrected::DepBool=nothing) +function var!(R::AbstractArray, A::AbstractArray{<:Real}, w::AbstractWeights, dims::Int; + mean=nothing, corrected::Union{Bool, Nothing}=nothing) corrected = depcheck(:var!, :corrected, corrected) if mean == 0 @@ -55,8 +55,8 @@ function var!(R::AbstractArray, A::RealArray, w::AbstractWeights, dims::Int; varcorrection(w, corrected)) end -function var(A::RealArray, w::AbstractWeights, dim::Int; mean=nothing, - corrected::DepBool=nothing) +function var(A::AbstractArray{<:Real}, w::AbstractWeights, dim::Int; mean=nothing, + corrected::Union{Bool, Nothing}=nothing) corrected = depcheck(:var, :corrected, corrected) var!(similar(A, Float64, Base.reduced_indices(axes(A), dim)), A, w, dim; mean=mean, corrected=corrected) @@ -81,11 +81,11 @@ weights used: * `ProbabilityWeights`: ``\\frac{n}{(n - 1) \\sum w}`` where ``n`` equals `count(!iszero, w)` * `Weights`: `ArgumentError` (bias correction not supported) """ -std(v::RealArray, w::AbstractWeights; mean=nothing, corrected::DepBool=nothing) = +std(v::AbstractArray{<:Real}, w::AbstractWeights; mean=nothing, corrected::Union{Bool, Nothing}=nothing) = sqrt.(var(v, w; mean=mean, corrected=depcheck(:std, :corrected, corrected))) -std(v::RealArray, w::AbstractWeights, dim::Int; - mean=nothing, corrected::DepBool=nothing) = +std(v::AbstractArray{<:Real}, w::AbstractWeights, dim::Int; + mean=nothing, corrected::Union{Bool, Nothing}=nothing) = sqrt.(var(v, w, dim; mean=mean, corrected=depcheck(:std, :corrected, corrected))) ##### Fused statistics @@ -120,38 +120,38 @@ function mean_and_std(x; corrected::Bool=true) m, s end -function mean_and_var(x::RealArray, w::AbstractWeights; corrected::DepBool=nothing) +function mean_and_var(x::AbstractArray{<:Real}, w::AbstractWeights; corrected::Union{Bool, Nothing}=nothing) m = mean(x, w) v = var(x, w, mean=m, corrected=depcheck(:mean_and_var, :corrected, corrected)) m, v end -function mean_and_std(x::RealArray, w::AbstractWeights; corrected::DepBool=nothing) +function mean_and_std(x::AbstractArray{<:Real}, w::AbstractWeights; corrected::Union{Bool, Nothing}=nothing) m = mean(x, w) s = std(x, w, mean=m, corrected=depcheck(:mean_and_std, :corrected, corrected)) m, s end -function mean_and_var(x::RealArray, dim::Int; corrected::Bool=true) +function mean_and_var(x::AbstractArray{<:Real}, dim::Int; corrected::Bool=true) m = mean(x, dims=dim) v = var(x, dims=dim, mean=m, corrected=corrected) m, v end -function mean_and_std(x::RealArray, dim::Int; corrected::Bool=true) +function mean_and_std(x::AbstractArray{<:Real}, dim::Int; corrected::Bool=true) m = mean(x, dims=dim) s = std(x, dims=dim, mean=m, corrected=corrected) m, s end -function mean_and_var(x::RealArray, w::AbstractWeights, dims::Int; - corrected::DepBool=nothing) +function mean_and_var(x::AbstractArray{<:Real}, w::AbstractWeights, dims::Int; + corrected::Union{Bool, Nothing}=nothing) m = mean(x, w, dims=dims) v = var(x, w, dims, mean=m, corrected=depcheck(:mean_and_var, :corrected, corrected)) m, v end -function mean_and_std(x::RealArray, w::AbstractWeights, dims::Int; - corrected::DepBool=nothing) +function mean_and_std(x::AbstractArray{<:Real}, w::AbstractWeights, dims::Int; + corrected::Union{Bool, Nothing}=nothing) m = mean(x, w, dims=dims) s = std(x, w, dims, mean=m, corrected=depcheck(:mean_and_std, :corrected, corrected)) m, s @@ -160,7 +160,7 @@ end ##### General central moment -function _moment2(v::RealArray, m::Real; corrected=false) +function _moment2(v::AbstractArray{<:Real}, m::Real; corrected=false) n = length(v) s = 0.0 for i = 1:n @@ -170,7 +170,7 @@ 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::AbstractArray{<:Real}, wv::AbstractWeights, m::Real; corrected=false) n = length(v) s = 0.0 for i = 1:n @@ -181,7 +181,7 @@ function _moment2(v::RealArray, wv::AbstractWeights, m::Real; corrected=false) varcorrection(wv, corrected) * s end -function _moment3(v::RealArray, m::Real) +function _moment3(v::AbstractArray{<:Real}, m::Real) n = length(v) s = 0.0 for i = 1:n @@ -191,7 +191,7 @@ function _moment3(v::RealArray, m::Real) s / n end -function _moment3(v::RealArray, wv::AbstractWeights, m::Real) +function _moment3(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) n = length(v) s = 0.0 for i = 1:n @@ -201,7 +201,7 @@ function _moment3(v::RealArray, wv::AbstractWeights, m::Real) s / sum(wv) end -function _moment4(v::RealArray, m::Real) +function _moment4(v::AbstractArray{<:Real}, m::Real) n = length(v) s = 0.0 for i = 1:n @@ -211,7 +211,7 @@ function _moment4(v::RealArray, m::Real) s / n end -function _moment4(v::RealArray, wv::AbstractWeights, m::Real) +function _moment4(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) n = length(v) s = 0.0 for i = 1:n @@ -221,7 +221,7 @@ function _moment4(v::RealArray, wv::AbstractWeights, m::Real) s / sum(wv) end -function _momentk(v::RealArray, k::Int, m::Real) +function _momentk(v::AbstractArray{<:Real}, k::Int, m::Real) n = length(v) s = 0.0 for i = 1:n @@ -231,7 +231,7 @@ function _momentk(v::RealArray, k::Int, m::Real) s / n end -function _momentk(v::RealArray, k::Int, wv::AbstractWeights, m::Real) +function _momentk(v::AbstractArray{<:Real}, k::Int, wv::AbstractWeights, m::Real) n = length(v) s = 0.0 for i = 1:n @@ -248,22 +248,22 @@ end Return the `k`th order central moment of a real-valued array `v`, optionally specifying a weighting vector `wv` and a center `m`. """ -function moment(v::RealArray, k::Int, m::Real) +function moment(v::AbstractArray{<:Real}, k::Int, m::Real) k == 2 ? _moment2(v, m) : k == 3 ? _moment3(v, m) : k == 4 ? _moment4(v, m) : _momentk(v, k, m) end -function moment(v::RealArray, k::Int, wv::AbstractWeights, m::Real) +function moment(v::AbstractArray{<:Real}, k::Int, wv::AbstractWeights, m::Real) k == 2 ? _moment2(v, wv, m) : k == 3 ? _moment3(v, wv, m) : k == 4 ? _moment4(v, wv, m) : _momentk(v, k, wv, m) end -moment(v::RealArray, k::Int) = moment(v, k, mean(v)) -function moment(v::RealArray, k::Int, wv::AbstractWeights) +moment(v::AbstractArray{<:Real}, k::Int) = moment(v, k, mean(v)) +function moment(v::AbstractArray{<:Real}, k::Int, wv::AbstractWeights) moment(v, k, wv, mean(v, wv)) end @@ -278,7 +278,7 @@ 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::AbstractArray{<:Real}, m::Real) n = length(v) cm2 = 0.0 # empirical 2nd centered moment (variance) cm3 = 0.0 # empirical 3rd centered moment @@ -294,7 +294,7 @@ 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::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) cm2 = 0.0 # empirical 2nd centered moment (variance) @@ -314,8 +314,8 @@ function skewness(v::RealArray, wv::AbstractWeights, m::Real) return cm3 / sqrt(cm2 * cm2 * cm2) # this is much faster than cm2^1.5 end -skewness(v::RealArray) = skewness(v, mean(v)) -skewness(v::RealArray, wv::AbstractWeights) = skewness(v, wv, mean(v, wv)) +skewness(v::AbstractArray{<:Real}) = skewness(v, mean(v)) +skewness(v::AbstractArray{<:Real}, wv::AbstractWeights) = skewness(v, wv, mean(v, wv)) # (excessive) Kurtosis # This is Type 1 definition according to Joanes and Gill (1998) @@ -325,7 +325,7 @@ 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::AbstractArray{<:Real}, m::Real) n = length(v) cm2 = 0.0 # empirical 2nd centered moment (variance) cm4 = 0.0 # empirical 4th centered moment @@ -340,7 +340,7 @@ function kurtosis(v::RealArray, m::Real) return (cm4 / (cm2 * cm2)) - 3.0 end -function kurtosis(v::RealArray, wv::AbstractWeights, m::Real) +function kurtosis(v::AbstractArray{<:Real}, wv::AbstractWeights, m::Real) n = length(v) length(wv) == n || throw(DimensionMismatch("Inconsistent array lengths.")) cm2 = 0.0 # empirical 2nd centered moment (variance) @@ -361,8 +361,8 @@ function kurtosis(v::RealArray, wv::AbstractWeights, m::Real) return (cm4 / (cm2 * cm2)) - 3.0 end -kurtosis(v::RealArray) = kurtosis(v, mean(v)) -kurtosis(v::RealArray, wv::AbstractWeights) = kurtosis(v, wv, mean(v, wv)) +kurtosis(v::AbstractArray{<:Real}) = kurtosis(v, mean(v)) +kurtosis(v::AbstractArray{<:Real}, wv::AbstractWeights) = kurtosis(v, wv, mean(v, wv)) """ cumulant(v, k, [wv::AbstractWeights], m=mean(v)) @@ -378,7 +378,7 @@ Reference: Smith, P. J. 1995. A Recursive Formulation of the Old Problem of Obta Moments from Cumulants and Vice Versa. The American Statistician, 49(2), 217–218. https://doi.org/10.2307/2684642 """ -function cumulant(v::RealArray, krange::Union{Integer, AbstractRange{<:Integer}}, wv::AbstractWeights, +function cumulant(v::AbstractArray{<:Real}, krange::Union{Integer, AbstractRange{<:Integer}}, wv::AbstractWeights, m::Real=mean(v, wv)) if minimum(krange) <= 0 throw(ArgumentError("Cumulant orders must be strictly positive.")) @@ -399,5 +399,5 @@ function cumulant(v::RealArray, krange::Union{Integer, AbstractRange{<:Integer}} return cumls[krange] end -cumulant(v::RealArray, krange::Union{Integer, AbstractRange{<:Integer}}, m::Real=mean(v)) = +cumulant(v::AbstractArray{<:Real}, krange::Union{Integer, AbstractRange{<:Integer}}, m::Real=mean(v)) = cumulant(v, krange, uweights(length(v)), m) diff --git a/src/rankcorr.jl b/src/rankcorr.jl index 00eee4380..4e023a018 100644 --- a/src/rankcorr.jl +++ b/src/rankcorr.jl @@ -17,14 +17,14 @@ Compute Spearman's rank correlation coefficient. If `x` and `y` are vectors, the output is a float, otherwise it's a matrix corresponding to the pairwise correlations of the columns of `x` and `y`. """ -function corspearman(x::RealVector, y::RealVector) +function corspearman(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) n = length(x) n == length(y) || throw(DimensionMismatch("vectors must have same length")) (any(isnan, x) || any(isnan, y)) && return NaN return cor(tiedrank(x), tiedrank(y)) end -function corspearman(X::RealMatrix, y::RealVector) +function corspearman(X::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}) size(X, 1) == length(y) || throw(DimensionMismatch("X and y have inconsistent dimensions")) n = size(X, 2) @@ -43,7 +43,7 @@ function corspearman(X::RealMatrix, y::RealVector) return C end -function corspearman(x::RealVector, Y::RealMatrix) +function corspearman(x::AbstractVector{<:Real}, Y::AbstractMatrix{<:Real}) size(Y, 1) == length(x) || throw(DimensionMismatch("x and Y have inconsistent dimensions")) n = size(Y, 2) @@ -62,7 +62,7 @@ function corspearman(x::RealVector, Y::RealMatrix) return C end -function corspearman(X::RealMatrix) +function corspearman(X::AbstractMatrix{<:Real}) n = size(X, 2) C = Matrix{Float64}(I, n, n) anynan = Vector{Bool}(undef, n) @@ -89,7 +89,7 @@ function corspearman(X::RealMatrix) return C end -function corspearman(X::RealMatrix, Y::RealMatrix) +function corspearman(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}) size(X, 1) == size(Y, 1) || throw(ArgumentError("number of rows in each array must match")) nr = size(X, 2) @@ -125,7 +125,7 @@ end # Knight, William R. “A Computer Method for Calculating Kendall's Tau with Ungrouped Data.” # Journal of the American Statistical Association, vol. 61, no. 314, 1966, pp. 436–439. # JSTOR, www.jstor.org/stable/2282833. -function corkendall!(x::RealVector, y::RealVector, permx::AbstractVector{<:Integer}=sortperm(x)) +function corkendall!(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, permx::AbstractArray{<:Integer}=sortperm(x)) if any(isnan, x) || any(isnan, y) return NaN end n = length(x) if n != length(y) error("Vectors must have same length") end @@ -173,20 +173,20 @@ end Compute Kendall's rank correlation coefficient, τ. `x` and `y` must both be either matrices or vectors. """ -corkendall(x::RealVector, y::RealVector) = corkendall!(copy(x), copy(y)) +corkendall(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}) = corkendall!(copy(x), copy(y)) -function corkendall(X::RealMatrix, y::RealVector) +function corkendall(X::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}) permy = sortperm(y) return([corkendall!(copy(y), X[:,i], permy) for i in 1:size(X, 2)]) end -function corkendall(x::RealVector, Y::RealMatrix) +function corkendall(x::AbstractVector{<:Real}, Y::AbstractMatrix{<:Real}) n = size(Y, 2) permx = sortperm(x) return(reshape([corkendall!(copy(x), Y[:,i], permx) for i in 1:n], 1, n)) end -function corkendall(X::RealMatrix) +function corkendall(X::AbstractMatrix{<:Real}) n = size(X, 2) C = Matrix{Float64}(I, n, n) for j = 2:n @@ -199,7 +199,7 @@ function corkendall(X::RealMatrix) return C end -function corkendall(X::RealMatrix, Y::RealMatrix) +function corkendall(X::AbstractMatrix{<:Real}, Y::AbstractMatrix{<:Real}) nr = size(X, 2) nc = size(Y, 2) C = Matrix{Float64}(undef, nr, nc) @@ -215,7 +215,7 @@ end # Auxilliary functions for Kendall's rank correlation """ - countties(x::RealVector, lo::Integer, hi::Integer) + countties(x::AbstractVector{<:Real}, lo::Integer, hi::Integer) Return the number of ties within `x[lo:hi]`. Assumes `x` is sorted. """ diff --git a/src/ranking.jl b/src/ranking.jl index 05a5b4657..cc0418894 100644 --- a/src/ranking.jl +++ b/src/ranking.jl @@ -31,7 +31,7 @@ function _rank(f!, x::AbstractArray{>: Missing}, R::Type=Int; sortkwargs...) end # Ordinal ranking ("1234 ranking") -- use the literal order resulted from sort -function _ordinalrank!(rks::AbstractArray, x::AbstractArray, p::IntegerArray) +function _ordinalrank!(rks::AbstractArray, x::AbstractArray, p::AbstractArray{<:Integer}) _check_randparams(rks, x, p) @inbounds for i in eachindex(p) rks[p[i]] = i @@ -54,7 +54,7 @@ ordinalrank(x::AbstractArray; sortkwargs...) = # Competition ranking ("1224" ranking) -- resolve tied ranks using min -function _competerank!(rks::AbstractArray, x::AbstractArray, p::IntegerArray) +function _competerank!(rks::AbstractArray, x::AbstractArray, p::AbstractArray{<:Integer}) n = _check_randparams(rks, x, p) @inbounds if n > 0 @@ -91,7 +91,7 @@ competerank(x::AbstractArray; sortkwargs...) = # Dense ranking ("1223" ranking) -- resolve tied ranks using min -function _denserank!(rks::AbstractArray, x::AbstractArray, p::IntegerArray) +function _denserank!(rks::AbstractArray, x::AbstractArray, p::AbstractArray{<:Integer}) n = _check_randparams(rks, x, p) @inbounds if n > 0 @@ -128,7 +128,7 @@ denserank(x::AbstractArray; sortkwargs...) = # Tied ranking ("1 2.5 2.5 4" ranking) -- resolve tied ranks using average -function _tiedrank!(rks::AbstractArray, x::AbstractArray, p::IntegerArray) +function _tiedrank!(rks::AbstractArray, x::AbstractArray, p::AbstractArray{<:Integer}) n = _check_randparams(rks, x, p) @inbounds if n > 0 diff --git a/src/sampling.jl b/src/sampling.jl index ea19a9306..0e5662e71 100644 --- a/src/sampling.jl +++ b/src/sampling.jl @@ -1028,10 +1028,10 @@ items appear in the same order as in `a`) should be taken. Optionally specify a random number generator `rng` as the first argument (defaults to `Random.GLOBAL_RNG`). """ -wsample!(rng::AbstractRNG, a::AbstractArray, w::RealVector, x::AbstractArray; +wsample!(rng::AbstractRNG, a::AbstractArray, w::AbstractVector{<:Real}, x::AbstractArray; replace::Bool=true, ordered::Bool=false) = sample!(rng, a, weights(w), x; replace=replace, ordered=ordered) -wsample!(a::AbstractArray, w::RealVector, x::AbstractArray; +wsample!(a::AbstractArray, w::AbstractVector{<:Real}, x::AbstractArray; replace::Bool=true, ordered::Bool=false) = sample!(Random.GLOBAL_RNG, a, weights(w), x; replace=replace, ordered=ordered) @@ -1044,10 +1044,10 @@ to the weights given in `w`. If `a` is not present, select a random weight from Optionally specify a random number generator `rng` as the first argument (defaults to `Random.GLOBAL_RNG`). """ -wsample(rng::AbstractRNG, w::RealVector) = sample(rng, weights(w)) -wsample(w::RealVector) = wsample(Random.GLOBAL_RNG, w) -wsample(rng::AbstractRNG, a::AbstractArray, w::RealVector) = sample(rng, a, weights(w)) -wsample(a::AbstractArray, w::RealVector) = wsample(Random.GLOBAL_RNG, a, w) +wsample(rng::AbstractRNG, w::AbstractVector{<:Real}) = sample(rng, weights(w)) +wsample(w::AbstractVector{<:Real}) = wsample(Random.GLOBAL_RNG, w) +wsample(rng::AbstractRNG, a::AbstractArray, w::AbstractVector{<:Real}) = sample(rng, a, weights(w)) +wsample(a::AbstractArray, w::AbstractVector{<:Real}) = wsample(Random.GLOBAL_RNG, a, w) """ @@ -1063,10 +1063,10 @@ items appear in the same order as in `a`) should be taken. Optionally specify a random number generator `rng` as the first argument (defaults to `Random.GLOBAL_RNG`). """ -wsample(rng::AbstractRNG, a::AbstractArray{T}, w::RealVector, n::Integer; +wsample(rng::AbstractRNG, a::AbstractArray{T}, w::AbstractVector{<:Real}, n::Integer; replace::Bool=true, ordered::Bool=false) where {T} = wsample!(rng, a, w, Vector{T}(undef, n); replace=replace, ordered=ordered) -wsample(a::AbstractArray, w::RealVector, n::Integer; +wsample(a::AbstractArray, w::AbstractVector{<:Real}, n::Integer; replace::Bool=true, ordered::Bool=false) = wsample(Random.GLOBAL_RNG, a, w, n; replace=replace, ordered=ordered) @@ -1080,9 +1080,9 @@ weights given in `w` if `a` is present, otherwise select a random sample of size Optionally specify a random number generator `rng` as the first argument (defaults to `Random.GLOBAL_RNG`). """ -wsample(rng::AbstractRNG, a::AbstractArray{T}, w::RealVector, dims::Dims; +wsample(rng::AbstractRNG, a::AbstractArray{T}, w::AbstractVector{<:Real}, dims::Dims; replace::Bool=true, ordered::Bool=false) where {T} = wsample!(rng, a, w, Array{T}(undef, dims); replace=replace, ordered=ordered) -wsample(a::AbstractArray, w::RealVector, dims::Dims; +wsample(a::AbstractArray, w::AbstractVector{<:Real}, dims::Dims; replace::Bool=true, ordered::Bool=false) = wsample(Random.GLOBAL_RNG, a, w, dims; replace=replace, ordered=ordered) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index bc3bfe0fd..6596c78a5 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -14,7 +14,7 @@ default_laglen(lx::Int) = min(lx-1, round(Int,10*log10(lx))) check_lags(lx::Int, lags::AbstractVector) = (maximum(lags) < lx || error("lags must be less than the sample length.")) -function demean_col!(z::RealVector, x::RealMatrix, j::Int, demean::Bool) +function demean_col!(z::AbstractVector{<:Real}, x::AbstractMatrix{<:Real}, j::Int, demean::Bool) T = eltype(z) m = size(x, 1) @assert m == length(z) @@ -43,8 +43,8 @@ end default_autolags(lx::Int) = 0 : default_laglen(lx) -_autodot(x::AbstractVector{<:RealFP}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx) -_autodot(x::RealVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx)) +_autodot(x::AbstractVector{<:Union{Float32, Float64}}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx) +_autodot(x::AbstractVector{<:Real}, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx)) ## autocov @@ -61,7 +61,7 @@ where each column in the result will correspond to a column in `x`. The output is not normalized. See [`autocor!`](@ref) for a method with normalization. """ -function autocov!(r::RealVector, x::RealVector, lags::IntegerVector; demean::Bool=true) +function autocov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) length(r) == m || throw(DimensionMismatch()) @@ -75,7 +75,7 @@ function autocov!(r::RealVector, x::RealVector, lags::IntegerVector; demean::Boo return r end -function autocov!(r::RealMatrix, x::RealMatrix, lags::IntegerVector; demean::Bool=true) +function autocov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -110,12 +110,12 @@ When left unspecified, the lags used are the integers from 0 to The output is not normalized. See [`autocor`](@ref) for a function with normalization. """ -function autocov(x::RealVector, lags::IntegerVector; demean::Bool=true) +function autocov(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) autocov!(out, x, lags; demean=demean) end -function autocov(x::RealMatrix, lags::IntegerVector; demean::Bool=true) +function autocov(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) autocov!(out, x, lags; demean=demean) end @@ -139,7 +139,7 @@ where each column in the result will correspond to a column in `x`. The output is normalized by the variance of `x`, i.e. so that the lag 0 autocorrelation is 1. See [`autocov!`](@ref) for the unnormalized form. """ -function autocor!(r::RealVector, x::RealVector, lags::IntegerVector; demean::Bool=true) +function autocor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) length(r) == m || throw(DimensionMismatch()) @@ -154,7 +154,7 @@ function autocor!(r::RealVector, x::RealVector, lags::IntegerVector; demean::Boo return r end -function autocor!(r::RealMatrix, x::RealMatrix, lags::IntegerVector; demean::Bool=true) +function autocor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -191,12 +191,12 @@ When left unspecified, the lags used are the integers from 0 to The output is normalized by the variance of `x`, i.e. so that the lag 0 autocorrelation is 1. See [`autocov`](@ref) for the unnormalized form. """ -function autocor(x::RealVector, lags::IntegerVector; demean::Bool=true) +function autocor(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) autocor!(out, x, lags; demean=demean) end -function autocor(x::RealMatrix, lags::IntegerVector; demean::Bool=true) +function autocor(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) autocor!(out, x, lags; demean=demean) end @@ -213,14 +213,14 @@ autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) = default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l) -function _crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:RealFP} +function _crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:Union{Float32, Float64}} if l >= 0 dot(x, 1:(lx-l), y, (1+l):lx) else dot(x, (1-l):lx, y, 1:(lx+l)) end end -function _crossdot(x::RealVector, y::RealVector, lx::Int, l::Int) +function _crossdot(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lx::Int, l::Int) if l >= 0 dot(view(x, 1:(lx-l)), view(y, (1+l):lx)) else @@ -246,7 +246,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`. The output is not normalized. See [`crosscor!`](@ref) for a function with normalization. """ -function crosscov!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) (length(y) == lx && length(r) == m) || throw(DimensionMismatch()) @@ -262,7 +262,7 @@ function crosscov!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVec return r end -function crosscov!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -282,7 +282,7 @@ function crosscov!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVec return r end -function crosscov!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) ns = size(y, 2) m = length(lags) @@ -302,7 +302,7 @@ function crosscov!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVec return r end -function crosscov!(r::AbstractArray{<:Real,3}, x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscov!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) @@ -356,22 +356,22 @@ When left unspecified, the lags used are the integers from The output is not normalized. See [`crosscor`](@ref) for a function with normalization. """ -function crosscov(x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscov(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags)) crosscov!(out, x, y, lags; demean=demean) end -function crosscov(x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscov(x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2)) crosscov!(out, x, y, lags; demean=demean) end -function crosscov(x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscov(x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2)) crosscov!(out, x, y, lags; demean=demean) end -function crosscov(x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscov(x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2)) crosscov!(out, x, y, lags; demean=demean) end @@ -397,7 +397,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`. The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov!`](@ref) for the unnormalized form. """ -function crosscor!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) (length(y) == lx && length(r) == m) || throw(DimensionMismatch()) @@ -414,7 +414,7 @@ function crosscor!(r::RealVector, x::RealVector, y::RealVector, lags::IntegerVec return r end -function crosscor!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -436,7 +436,7 @@ function crosscor!(r::RealMatrix, x::RealMatrix, y::RealVector, lags::IntegerVec return r end -function crosscor!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) ns = size(y, 2) m = length(lags) @@ -458,7 +458,7 @@ function crosscor!(r::RealMatrix, x::RealVector, y::RealMatrix, lags::IntegerVec return r end -function crosscor!(r::AbstractArray{<:Real,3}, x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscor!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) @@ -517,22 +517,22 @@ When left unspecified, the lags used are the integers from The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov`](@ref) for the unnormalized form. """ -function crosscor(x::RealVector, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscor(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags)) crosscor!(out, x, y, lags; demean=demean) end -function crosscor(x::RealMatrix, y::RealVector, lags::IntegerVector; demean::Bool=true) +function crosscor(x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2)) crosscor!(out, x, y, lags; demean=demean) end -function crosscor(x::RealVector, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscor(x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2)) crosscor!(out, x, y, lags; demean=demean) end -function crosscor(x::RealMatrix, y::RealMatrix, lags::IntegerVector; demean::Bool=true) +function crosscor(x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2)) crosscor!(out, x, y, lags; demean=demean) end @@ -549,7 +549,7 @@ crosscor(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool= # ####################################### -function pacf_regress!(r::RealMatrix, X::RealMatrix, lags::IntegerVector, mk::Integer) +function pacf_regress!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}, mk::Integer) lx = size(X, 1) tmpX = ones(eltype(X), lx, mk + 1) for j = 1 : size(X,2) @@ -561,13 +561,13 @@ function pacf_regress!(r::RealMatrix, X::RealMatrix, lags::IntegerVector, mk::In for i = 1 : length(lags) l = lags[i] sX = view(tmpX, 1+l:lx, 1:l+1) - r[i,j] = l == 0 ? 1 : (cholesky!(sX'sX, Val(false)) \ (sX'view(X, 1+l:lx, j)))[end] + r[i,j] = l == 0 ? 1 : (cholesky!(sX'sX) \ (sX'view(X, 1+l:lx, j)))[end] end end r end -function pacf_yulewalker!(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector, mk::Integer) where T<:RealFP +function pacf_yulewalker!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}, mk::Integer) where T<:Union{Float32, Float64} tmp = Vector{T}(undef, mk) for j = 1 : size(X,2) acfs = autocor(X[:,j], 1:mk) @@ -590,7 +590,7 @@ using the Yule-Walker equations. `r` must be a matrix of size `(length(lags), size(x, 2))`. """ -function pacf!(r::RealMatrix, X::AbstractMatrix{T}, lags::IntegerVector; method::Symbol=:regression) where T<:RealFP +function pacf!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) where T<:Union{Float32, Float64} lx = size(X, 1) m = length(lags) minlag, maxlag = extrema(lags) @@ -621,11 +621,11 @@ If `x` is a vector, return a vector of the same length as `lags`. If `x` is a matrix, return a matrix of size `(length(lags), size(x, 2))`, where each column in the result corresponds to a column in `x`. """ -function pacf(X::RealMatrix, lags::IntegerVector; method::Symbol=:regression) +function pacf(X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) out = Matrix{float(eltype(X))}(undef, length(lags), size(X,2)) pacf!(out, float(X), lags; method=method) end -function pacf(x::RealVector, lags::IntegerVector; method::Symbol=:regression) +function pacf(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) vec(pacf(reshape(x, length(x), 1), lags, method=method)) end diff --git a/src/weights.jl b/src/weights.jl index 78091c2ae..cf535d408 100644 --- a/src/weights.jl +++ b/src/weights.jl @@ -5,7 +5,7 @@ abstract type AbstractWeights{S<:Real, T<:Real, V<:AbstractVector{T}} <: Abstrac @weights name Generates a new generic weight type with specified `name`, which subtypes `AbstractWeights` -and stores the `values` (`V<:RealVector`) and `sum` (`S<:Real`). +and stores the `values` (`V<:AbstractVector{<:Real}`) and `sum` (`S<:Real`). """ macro weights(name) return quote @@ -85,8 +85,8 @@ and [`ProbabilityWeights`](@ref). Construct a `Weights` vector from array `vs`. See the documentation for [`Weights`](@ref) for more details. """ -weights(vs::RealArray) = Weights(vec(vs)) -weights(vs::RealVector) = Weights(vs) +weights(vs::AbstractArray{<:Real}) = Weights(vec(vs)) +weights(vs::AbstractVector{<:Real}) = Weights(vs) """ varcorrection(w::Weights, corrected=false) @@ -120,8 +120,8 @@ being weighted are aggregate values (e.g., averages) with differing variances. Construct an `AnalyticWeights` vector from array `vs`. See the documentation for [`AnalyticWeights`](@ref) for more details. """ -aweights(vs::RealVector) = AnalyticWeights(vs) -aweights(vs::RealArray) = AnalyticWeights(vec(vs)) +aweights(vs::AbstractVector{<:Real}) = AnalyticWeights(vs) +aweights(vs::AbstractArray{<:Real}) = AnalyticWeights(vec(vs)) """ varcorrection(w::AnalyticWeights, corrected=false) @@ -158,8 +158,8 @@ was observed. These weights may also be referred to as case weights or repeat we Construct a `FrequencyWeights` vector from a given array. See the documentation for [`FrequencyWeights`](@ref) for more details. """ -fweights(vs::RealVector) = FrequencyWeights(vs) -fweights(vs::RealArray) = FrequencyWeights(vec(vs)) +fweights(vs::AbstractVector{<:Real}) = FrequencyWeights(vs) +fweights(vs::AbstractArray{<:Real}) = FrequencyWeights(vec(vs)) """ varcorrection(w::FrequencyWeights, corrected=false) @@ -196,8 +196,8 @@ These weights may also be referred to as sampling weights. Construct a `ProbabilityWeights` vector from a given array. See the documentation for [`ProbabilityWeights`](@ref) for more details. """ -pweights(vs::RealVector) = ProbabilityWeights(vs) -pweights(vs::RealArray) = ProbabilityWeights(vec(vs)) +pweights(vs::AbstractVector{<:Real}) = ProbabilityWeights(vs) +pweights(vs::AbstractArray{<:Real}) = ProbabilityWeights(vec(vs)) """ varcorrection(w::ProbabilityWeights, corrected=false) @@ -217,7 +217,7 @@ pweights(vs::RealArray) = ProbabilityWeights(vec(vs)) end """ - eweights(t::AbstractVector{<:Integer}, λ::Real; scale=false) + eweights(t::AbstractArray{<:Integer}, λ::Real; scale=false) eweights(t::AbstractVector{T}, r::StepRange{T}, λ::Real; scale=false) where T eweights(n::Integer, λ::Real; scale=false) @@ -268,7 +268,7 @@ julia> eweights(1:10, 0.3; scale=true) - https://en.wikipedia.org/wiki/Moving_average#Exponential_moving_average - https://en.wikipedia.org/wiki/Exponential_smoothing """ -function eweights(t::AbstractVector{<:Integer}, λ::Real; kwargs...) +function eweights(t::AbstractArray{<:Integer}, λ::Real; kwargs...) isempty(t) && return Weights(copy(t), 0) (lo, hi) = extrema(t) return _eweights(t, λ, hi - lo + 1; kwargs...) @@ -278,7 +278,7 @@ eweights(n::Integer, λ::Real; kwargs...) = _eweights(1:n, λ, n; kwargs...) eweights(t::AbstractVector, r::AbstractRange, λ::Real; kwargs...) = _eweights(something.(indexin(t, r)), λ, length(r); kwargs...) -function _eweights(t::AbstractVector{<:Integer}, λ::Real, n::Integer; scale::DepBool=nothing) +function _eweights(t::AbstractArray{<:Integer}, λ::Real, n::Integer; scale::Union{Bool, Nothing}=nothing) 0 < λ <= 1 || throw(ArgumentError("Smoothing factor must be between 0 and 1")) f = depcheck(:eweights, :scale, scale) ? _scaled_eweight : _unscaled_eweight @@ -704,7 +704,7 @@ is strictly superior to ``h``. The weighted ``p`` quantile is given by ``v_k + \ with ``\\gamma = (h - S_k)/(S_{k+1} - S_k)``. In particular, when all weights are equal, the function returns the same result as the unweighted `quantile`. """ -function quantile(v::RealVector{V}, w::AbstractWeights{W}, p::RealVector) where {V,W<:Real} +function quantile(v::AbstractVector{<:Real}{V}, w::AbstractWeights{W}, p::AbstractVector{<:Real}) where {V,W<:Real} # checks isempty(v) && throw(ArgumentError("quantile of an empty array is undefined")) isempty(p) && throw(ArgumentError("empty quantile array")) @@ -770,19 +770,19 @@ function quantile(v::RealVector{V}, w::AbstractWeights{W}, p::RealVector) where return out end -function quantile(v::RealVector, w::UnitWeights, p::RealVector) +function quantile(v::AbstractVector{<:Real}, w::UnitWeights, p::AbstractVector{<:Real}) length(v) != length(w) && throw(DimensionMismatch("Inconsistent array dimension.")) return quantile(v, p) end -quantile(v::RealVector, w::AbstractWeights{<:Real}, p::Number) = quantile(v, w, [p])[1] +quantile(v::AbstractVector{<:Real}, w::AbstractWeights{<:Real}, p::Number) = quantile(v, w, [p])[1] ##### Weighted median ##### """ - median(v::RealVector, w::AbstractWeights) + median(v::AbstractVector{<:Real}, w::AbstractWeights) Compute the weighted median of `v` with weights `w` (of type `AbstractWeights`). See the documentation for [`quantile`](@ref) for more details. """ -median(v::RealVector, w::AbstractWeights{<:Real}) = quantile(v, w, 0.5) +median(v::AbstractVector{<:Real}, w::AbstractWeights{<:Real}) = quantile(v, w, 0.5) diff --git a/test/rankcorr.jl b/test/rankcorr.jl index 1e987a0ee..d1f58ef14 100644 --- a/test/rankcorr.jl +++ b/test/rankcorr.jl @@ -35,13 +35,13 @@ c22 = corspearman(x2, x2) @test isnan(corkendall([1,1,1], [1,2,3])) @test corkendall([-Inf,-0.0,Inf],[1,2,3]) == 1.0 -# Test, with exact equality, some known results. -# RealVector, RealVector +# Test, with exact equality, some known results. +# AbstractVector{<:Real}, AbstractVector{<:Real} @test corkendall(x1, y) == -1/sqrt(90) @test corkendall(x2, y) == -1/sqrt(72) -# RealMatrix, RealVector +# AbstractMatrix{<:Real}, AbstractVector{<:Real} @test corkendall(X, y) == [-1/sqrt(90), -1/sqrt(72)] -# RealVector, RealMatrix +# AbstractVector{<:Real}, AbstractMatrix{<:Real} @test corkendall(y, X) == [-1/sqrt(90) -1/sqrt(72)] # n = 78_000 tests for overflow errors on 32 bit @@ -66,9 +66,9 @@ c11 = corkendall(x1, x1) c12 = corkendall(x1, x2) c22 = corkendall(x2, x2) -# RealMatrix, RealMatrix +# AbstractMatrix{<:Real}, AbstractMatrix{<:Real} @test corkendall(X, X) ≈ [c11 c12; c12 c22] -# RealMatrix +# AbstractMatrix{<:Real} @test corkendall(X) ≈ [c11 c12; c12 c22] @test c11 == 1.0 @@ -158,4 +158,4 @@ end @test_throws ErrorException corkendall([1], [1, 2]) @test_throws ErrorException corkendall([1], [1 2; 3 4]) @test_throws ErrorException corkendall([1 2; 3 4], [1]) -@test_throws ArgumentError corkendall([1 2; 3 4: 4 6], [1 2; 3 4]) \ No newline at end of file +@test_throws ArgumentError corkendall([1 2; 3 4: 4 6], [1 2; 3 4]) From 3747453b627560013f3cbcc615aa279b8ca63602 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Sat, 27 Nov 2021 20:08:48 +0300 Subject: [PATCH 02/15] Compute auto- and cross-correlations for complex sequences --- src/signalcorr.jl | 69 +++++++++++---------- test/rcall_test.jl | 29 +++++++++ test/signalcorr.jl | 145 +++++++++++++++++++++++---------------------- 3 files changed, 136 insertions(+), 107 deletions(-) create mode 100644 test/rcall_test.jl diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 6596c78a5..3e6bf9e85 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -14,7 +14,7 @@ default_laglen(lx::Int) = min(lx-1, round(Int,10*log10(lx))) check_lags(lx::Int, lags::AbstractVector) = (maximum(lags) < lx || error("lags must be less than the sample length.")) -function demean_col!(z::AbstractVector{<:Real}, x::AbstractMatrix{<:Real}, j::Int, demean::Bool) +function demean_col!(z::AbstractVector, x::AbstractMatrix, j::Int, demean::Bool) T = eltype(z) m = size(x, 1) @assert m == length(z) @@ -43,8 +43,7 @@ end default_autolags(lx::Int) = 0 : default_laglen(lx) -_autodot(x::AbstractVector{<:Union{Float32, Float64}}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx) -_autodot(x::AbstractVector{<:Real}, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx)) +_autodot(x::AbstractVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx)) ## autocov @@ -61,7 +60,7 @@ where each column in the result will correspond to a column in `x`. The output is not normalized. See [`autocor!`](@ref) for a method with normalization. """ -function autocov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocov!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) length(r) == m || throw(DimensionMismatch()) @@ -75,7 +74,7 @@ function autocov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::Ab return r end -function autocov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocov!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -110,17 +109,17 @@ When left unspecified, the lags used are the integers from 0 to The output is not normalized. See [`autocor`](@ref) for a function with normalization. """ -function autocov(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocov(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) autocov!(out, x, lags; demean=demean) end -function autocov(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocov(x::AbstractMatrix, lags::AbstractVector; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) autocov!(out, x, lags; demean=demean) end -autocov(x::AbstractVecOrMat{<:Real}; demean::Bool=true) = +autocov(x::AbstractVecOrMat; demean::Bool=true) = autocov(x, default_autolags(size(x,1)); demean=demean) ## autocor @@ -139,7 +138,7 @@ where each column in the result will correspond to a column in `x`. The output is normalized by the variance of `x`, i.e. so that the lag 0 autocorrelation is 1. See [`autocov!`](@ref) for the unnormalized form. """ -function autocor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocor!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) length(r) == m || throw(DimensionMismatch()) @@ -154,7 +153,7 @@ function autocor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, lags::Ab return r end -function autocor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocor!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -191,17 +190,17 @@ When left unspecified, the lags used are the integers from 0 to The output is normalized by the variance of `x`, i.e. so that the lag 0 autocorrelation is 1. See [`autocov`](@ref) for the unnormalized form. """ -function autocor(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocor(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) autocor!(out, x, lags; demean=demean) end -function autocor(x::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocor(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) autocor!(out, x, lags; demean=demean) end -autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) = +autocor(x::AbstractVecOrMat; demean::Bool=true) = autocor(x, default_autolags(size(x,1)); demean=demean) @@ -213,14 +212,14 @@ autocor(x::AbstractVecOrMat{<:Real}; demean::Bool=true) = default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l) -function _crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:Union{Float32, Float64}} +function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int) if l >= 0 dot(x, 1:(lx-l), y, (1+l):lx) else dot(x, (1-l):lx, y, 1:(lx+l)) end end -function _crossdot(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lx::Int, l::Int) +function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int) if l >= 0 dot(view(x, 1:(lx-l)), view(y, (1+l):lx)) else @@ -246,7 +245,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`. The output is not normalized. See [`crosscor!`](@ref) for a function with normalization. """ -function crosscov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov!(r::AbstractVector, x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) (length(y) == lx && length(r) == m) || throw(DimensionMismatch()) @@ -262,7 +261,7 @@ function crosscov!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::Abst return r end -function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov!(r::AbstractMatrix, x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -282,7 +281,7 @@ function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::Abst return r end -function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) ns = size(y, 2) m = length(lags) @@ -302,7 +301,7 @@ function crosscov!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::Abst return r end -function crosscov!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov!(r::AbstractArray{U,3} where U, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) @@ -356,27 +355,27 @@ When left unspecified, the lags used are the integers from The output is not normalized. See [`crosscor`](@ref) for a function with normalization. """ -function crosscov(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov(x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags)) crosscov!(out, x, y, lags; demean=demean) end -function crosscov(x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov(x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2)) crosscov!(out, x, y, lags; demean=demean) end -function crosscov(x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov(x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2)) crosscov!(out, x, y, lags; demean=demean) end -function crosscov(x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov(x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2)) crosscov!(out, x, y, lags; demean=demean) end -crosscov(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=true) = +crosscov(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) = crosscov(x, y, default_crosslags(size(x,1)); demean=demean) @@ -397,7 +396,7 @@ three-dimensional array of size `(length(lags), size(x, 2), size(y, 2))`. The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov!`](@ref) for the unnormalized form. """ -function crosscor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor!(r::AbstractVector, x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) (length(y) == lx && length(r) == m) || throw(DimensionMismatch()) @@ -414,7 +413,7 @@ function crosscor!(r::AbstractVector{<:Real}, x::AbstractVector{<:Real}, y::Abst return r end -function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor!(r::AbstractMatrix, x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) ns = size(x, 2) m = length(lags) @@ -436,7 +435,7 @@ function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractMatrix{<:Real}, y::Abst return r end -function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) ns = size(y, 2) m = length(lags) @@ -458,7 +457,7 @@ function crosscor!(r::AbstractMatrix{<:Real}, x::AbstractVector{<:Real}, y::Abst return r end -function crosscor!(r::AbstractArray{<:Real,3}, x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor!(r::AbstractArray{U,3} where U, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) @@ -517,27 +516,27 @@ When left unspecified, the lags used are the integers from The output is normalized by `sqrt(var(x)*var(y))`. See [`crosscov`](@ref) for the unnormalized form. """ -function crosscor(x::AbstractVector{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor(x::AbstractVector, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(Base.promote_eltype(x, y))}(undef, length(lags)) crosscor!(out, x, y, lags; demean=demean) end -function crosscor(x::AbstractMatrix{<:Real}, y::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor(x::AbstractMatrix, y::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(x,2)) crosscor!(out, x, y, lags; demean=demean) end -function crosscor(x::AbstractVector{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor(x::AbstractVector, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(Base.promote_eltype(x, y))}(undef, length(lags), size(y,2)) crosscor!(out, x, y, lags; demean=demean) end -function crosscor(x::AbstractMatrix{<:Real}, y::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor(x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Array{float(Base.promote_eltype(x, y)),3}(undef, length(lags), size(x,2), size(y,2)) crosscor!(out, x, y, lags; demean=demean) end -crosscor(x::AbstractVecOrMat{<:Real}, y::AbstractVecOrMat{<:Real}; demean::Bool=true) = +crosscor(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) = crosscor(x, y, default_crosslags(size(x,1)); demean=demean) @@ -621,11 +620,11 @@ If `x` is a vector, return a vector of the same length as `lags`. If `x` is a matrix, return a matrix of size `(length(lags), size(x, 2))`, where each column in the result corresponds to a column in `x`. """ -function pacf(X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) +function pacf(X::RealMatrix, lags::IntegerVector; method::Symbol=:regression) out = Matrix{float(eltype(X))}(undef, length(lags), size(X,2)) pacf!(out, float(X), lags; method=method) end -function pacf(x::AbstractVector{<:Real}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) +function pacf(x::RealVector, lags::IntegerVector; method::Symbol=:regression) vec(pacf(reshape(x, length(x), 1), lags, method=method)) end diff --git a/test/rcall_test.jl b/test/rcall_test.jl new file mode 100644 index 000000000..efa0c1bfb --- /dev/null +++ b/test/rcall_test.jl @@ -0,0 +1,29 @@ +using RCall + +_r_ccf(x, y, lm) = rcopy(rcall(:ccf, x, y, plot=false, lag=lm, type="covariance"))[:acf] +_r_acf(x, lm) = rcopy(rcall(:acf, x, plot=false, lag=lm, type="covariance"))[:acf] +_r_pacf(x, lm) = rcopy(rcall(:acf, x, plot=false, lag=lm))[:acf] + +function _ccf(is_auto::Bool, x::AbstractVector{<:Real}...) + lx = size(x[1], 1) + lag_max = min(lx - 1, round(Int, 10*log10(lx))) + cor = reverse(dropdims(_r_ccf(x[1], x[end], lag_max), dims = (2, 3))) + is_auto ? cor[lx:end] : cor +end + +r_autocov(x::AbstractVector) = r_crosscov(x, x, true) + +function r_crosscov(x::AbstractVector, y::AbstractVector, is_auto=false) + cv(x,y) = _ccf(is_auto, x, y) + cc_re = cv(real.(x), real.(y)) + cv(imag.(x), imag.(y)) + cc_im = cv(real.(x), imag.(y)) - cv(imag.(x), real.(y)) + cc_re + im*cc_im +end + +r_autocor(x::AbstractVector) = r_crosscor(x, x, true) + +function r_crosscor(x::AbstractVector, y::AbstractVector, is_auto=false) + cc = r_crosscov(x, y, is_auto) + sc(z) = _r_ccf(real.(z), real.(z), 0) + _r_ccf(imag.(z), imag.(z), 0) + cc/sqrt(sc(x)*sc(y)) +end diff --git a/test/signalcorr.jl b/test/signalcorr.jl index bfbe90fed..3234ce060 100644 --- a/test/signalcorr.jl +++ b/test/signalcorr.jl @@ -9,88 +9,76 @@ using Test # random data for testing -x = [-2.133252557240862 -.7445937365828654; - .1775816414485478 -.5834801838041446; - -.6264517920318317 -.68444205333293; - -.8809042583216906 .9071671734302398; - .09251017186697393 -1.0404476733379926; - -.9271887119115569 -.620728578941385; - 3.355819743178915 -.8325051361909978; - -.2834039258495755 -.22394811874731657; - .5354280026977677 .7481337671592626; - .39182285417742585 .3085762550821047] +x = [0.5316732188954801 + 0.015032026010878084im -0.8051436916816284 - 0.730238071058262im + -0.26236085703683865 + 0.1561582684776427im -0.7733975759329091 + 0.7517417038318899im + -0.12750317530769195 + 0.8359780109572358im -0.1426497126379726 - 0.1976368614625515im + -0.3779012097765515 - 0.15588720329528907im 0.4625694946681199 - 0.6687893051682382im + -0.14198553599197566 - 0.911613885474317im -0.38003834312962825 - 0.16244828932805647im + 0.6809910683979451 - 1.1951424604960583im 0.36127363787490463 + 1.3578919531660747im + -0.7024909626802256 - 0.2915841698461933im 0.8432078641246675 - 1.373657703019187im + -0.5684618031428246 + 0.3729428437838556im -0.18853209977543475 - 0.9698921354142658im + 1.2992421672413557 - 0.17943637636017656im -1.7801079183259816 - 1.3063359194617465im + -0.26481424468298875 - 1.1636864449113618im -0.2605299632484854 - 0.7356290016171836im] x1 = view(x, :, 1) x2 = view(x, :, 2) -realx = convert(AbstractMatrix{Real}, x) -realx1 = convert(AbstractVector{Real}, x1) -realx2 = convert(AbstractVector{Real}, x2) # autocov & autocorr @test autocov([1:5;]) ≈ [2.0, 0.8, -0.2, -0.8, -0.8] @test autocor([1, 2, 3, 4, 5]) ≈ [1.0, 0.4, -0.1, -0.4, -0.4] -racovx1 = [1.839214242630635709475, - -0.406784553146903871124, - 0.421772254824993531042, - 0.035874943792884653182, - -0.255679775928512320604, - 0.231154400105831353551, - -0.787016960267425180753, - 0.039909287349160660341, - -0.110149697877911914579, - -0.088687020167434751916] +racovx1 = [0.755284179631112 + 0.0im, + -0.005333112170365584 - 0.18633291805458002im, + -0.28638133506011604 + 0.1397225175604478im, + -0.05476759020476188 + 0.12991227531087618im, + -0.00499902418309206 + 0.023463421186162473im, + 0.11854989361702409 - 0.0028772373657594066im, + -0.005123812003979093 - 0.08041878599400383im, + -0.14090692583679154 + 0.035230042290069444im, + 0.039899180711674635 + 0.004918236900383659im, + -0.03857936468514856 - 0.04063999068714769im] @test autocov(x1) ≈ racovx1 -@test autocov(realx1) ≈ racovx1 @test autocov(x) ≈ [autocov(x1) autocov(x2)] -@test autocov(realx) ≈ [autocov(realx1) autocov(realx2)] - -racorx1 = [0.999999999999999888978, - -0.221173011668873431557, - 0.229321981664153962122, - 0.019505581764945757045, - -0.139015765538446717242, - 0.125681062460244019618, - -0.427909344123907742219, - 0.021699096507690283225, - -0.059889541590524189574, - -0.048220059475281865091] + +racorx1 = [1.0 + 0.0im, + -0.007061066965510023 - 0.24670570770539213im, + -0.3791703080554228 + 0.18499330626611243im, + -0.07251256107537019 + 0.17200449686941222im, + -0.006618732813301651 + 0.031065686027770673im, + 0.1569606471499575 - 0.0038094765432061303im, + -0.00678395250709688 - 0.106474871528861im, + -0.18656146869859214 + 0.04664475073114351im, + 0.05282671316002114 + 0.0065117700502952056im, + -0.051079270194684986 - 0.0538075492419246im] @test autocor(x1) ≈ racorx1 -@test autocor(realx1) ≈ racorx1 @test autocor(x) ≈ [autocor(x1) autocor(x2)] -@test autocor(realx) ≈ [autocor(realx1) autocor(realx2)] - # crosscov & crosscor -rcov0 = [0.320000000000000006661, - -0.319999999999999951150, - 0.080000000000000029421, - -0.479999999999999982236, - 0.000000000000000000000, - 0.479999999999999982236, - -0.080000000000000029421, - 0.319999999999999951150, - -0.320000000000000006661] +rcov0 = [0.96 + 0.32im, + -0.96 - 0.32im, + 0.24 + 0.08im, + -1.44 - 0.48im, + 0.0 + 0.0im, + 1.44 + 0.48im, + -0.24 - 0.08im, + 0.96 + 0.32im, + -0.96 - 0.32im] -@test crosscov([1, 2, 3, 4, 5], [1, -1, 1, -1, 1]) ≈ rcov0 -@test crosscov([1:5;], [1:5;]) ≈ [-0.8, -0.8, -0.2, 0.8, 2.0, 0.8, -0.2, -0.8, -0.8] +@test crosscov([1+9im, 2+7im, 3+5im, 4+3im, 5+1im], [1-im, -1+im, 1-im, -1+im, 1-im]) ≈ rcov0 +@test crosscov([1:5 ;], [1:5;]) ≈ [-0.8, -0.8, -0.2, 0.8, 2.0, 0.8, -0.2, -0.8, -0.8] c11 = crosscov(x1, x1) c12 = crosscov(x1, x2) c21 = crosscov(x2, x1) c22 = crosscov(x2, x2) -@test crosscov(realx1, realx2) ≈ c12 @test crosscov(x, x1) ≈ [c11 c21] -@test crosscov(realx, realx1) ≈ [c11 c21] @test crosscov(x1, x) ≈ [c11 c12] -@test crosscov(realx1, realx) ≈ [c11 c12] @test crosscov(x, x) ≈ cat([c11 c21], [c12 c22], dims=3) -@test crosscov(realx, realx) ≈ cat([c11 c21], [c12 c22], dims=3) # issue #805: avoid converting one input to the other's eltype @test crosscov([34566.5345, 3466.4566], Float16[1, 10]) ≈ @@ -114,14 +102,10 @@ c11 = crosscor(x1, x1) c12 = crosscor(x1, x2) c21 = crosscor(x2, x1) c22 = crosscor(x2, x2) -@test crosscor(realx1, realx2) ≈ c12 @test crosscor(x, x1) ≈ [c11 c21] -@test crosscor(realx, realx1) ≈ [c11 c21] @test crosscor(x1, x) ≈ [c11 c12] -@test crosscor(realx1, realx) ≈ [c11 c12] @test crosscor(x, x) ≈ cat([c11 c21], [c12 c22], dims=3) -@test crosscor(realx, realx) ≈ cat([c11 c21], [c12 c22], dims=3) # issue #805: avoid converting one input to the other's eltype @test crosscor([34566.5345, 3466.4566], Float16[1, 10]) ≈ @@ -130,17 +114,34 @@ c22 = crosscor(x2, x2) ## pacf - -rpacfr = [-0.218158122381419, - 0.195015316828711, - 0.144315804606139, - -0.199791229449779] - -@test pacf(x[:,1], 1:4) ≈ rpacfr - -rpacfy = [-0.221173011668873, - 0.189683314308021, - 0.111857020733719, - -0.175020669835420] - -@test pacf(x[:,1], 1:4, method=:yulewalker) ≈ rpacfy +x_pacf = [-2.133252557240862 + .1775816414485478 + -.6264517920318317 + -.8809042583216906 + .09251017186697393 + -.9271887119115569 + 3.355819743178915 + -.2834039258495755 + .5354280026977677 + .39182285417742585] + +rpacfr = [-0.2181581223814188 + 0.1950153168287112 + 0.1443158046061393 + -0.1997912294497786 + 1.686452299291934 + -1.446624527965185] + +@test pacf(real.(x_pacf), 1:4) ≈ rpacfr[1:4] + +rpacfy = [-0.2211730116688735 + 0.1896833143080208 + 0.1118570207337193 + -0.1750206698354201 + 0.04490978543108597 + -0.3715806391613901 + -0.1743836518336153 + 0.07574637899230210 + 0.04970724340234131] + +@test pacf(real.(x_pacf), 1:4, method=:yulewalker) ≈ rpacfy[1:4] From b62800cbf09a6e99367fe537121b3958bf2fb770 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Mon, 29 Nov 2021 19:49:56 +0300 Subject: [PATCH 03/15] Durbin algorithm works correctly --- src/signalcorr.jl | 13 ++++++------ src/toeplitzsolvers.jl | 12 +++++------ test/signalcorr.jl | 45 +++++++++++++----------------------------- 3 files changed, 26 insertions(+), 44 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 3e6bf9e85..51e4cb1d0 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -65,7 +65,6 @@ function autocov!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:I m = length(lags) length(r) == m || throw(DimensionMismatch()) check_lags(lx, lags) - T = typeof(zero(eltype(x)) / 1) z::Vector{T} = demean ? x .- mean(x) : x for k = 1 : m # foreach lag value @@ -160,7 +159,7 @@ function autocor!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:I size(r) == (m, ns) || throw(DimensionMismatch()) check_lags(lx, lags) - T = typeof(zero(eltype(x)) / 1) + T = typeof(zero(eltype(x))/1) z = Vector{T}(undef, lx) for j = 1 : ns demean_col!(z, x, j, demean) @@ -548,7 +547,7 @@ crosscor(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) = # ####################################### -function pacf_regress!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}, lags::AbstractVector{<:Integer}, mk::Integer) +function pacf_regress!(r::AbstractMatrix X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer) lx = size(X, 1) tmpX = ones(eltype(X), lx, mk + 1) for j = 1 : size(X,2) @@ -566,7 +565,7 @@ function pacf_regress!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{<:Real}, lag r end -function pacf_yulewalker!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}, mk::Integer) where T<:Union{Float32, Float64} +function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer) where T<:Union{Float32, Float64} tmp = Vector{T}(undef, mk) for j = 1 : size(X,2) acfs = autocor(X[:,j], 1:mk) @@ -589,7 +588,7 @@ using the Yule-Walker equations. `r` must be a matrix of size `(length(lags), size(x, 2))`. """ -function pacf!(r::AbstractMatrix{<:Real}, X::AbstractMatrix{T}, lags::AbstractVector{<:Integer}; method::Symbol=:regression) where T<:Union{Float32, Float64} +function pacf!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector; method::Symbol=:regression) lx = size(X, 1) m = length(lags) minlag, maxlag = extrema(lags) @@ -620,11 +619,11 @@ If `x` is a vector, return a vector of the same length as `lags`. If `x` is a matrix, return a matrix of size `(length(lags), size(x, 2))`, where each column in the result corresponds to a column in `x`. """ -function pacf(X::RealMatrix, lags::IntegerVector; method::Symbol=:regression) +function pacf(X::AbstractMatrix, lags::IntegerVector; method::Symbol=:regression) out = Matrix{float(eltype(X))}(undef, length(lags), size(X,2)) pacf!(out, float(X), lags; method=method) end -function pacf(x::RealVector, lags::IntegerVector; method::Symbol=:regression) +function pacf(x::AbstractVector, lags::IntegerVector; method::Symbol=:regression) vec(pacf(reshape(x, length(x), 1), lags, method=method)) end diff --git a/src/toeplitzsolvers.jl b/src/toeplitzsolvers.jl index 19146bf7a..28ff92423 100644 --- a/src/toeplitzsolvers.jl +++ b/src/toeplitzsolvers.jl @@ -1,12 +1,12 @@ # Symmetric Toeplitz solver -function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T<:BlasReal +function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T n = length(r) n <= length(y) || throw(DimensionMismatch("Auxiliary vector cannot be shorter than data vector")) y[1] = -r[1] β = one(T) α = -r[1] for k = 1:n-1 - β *= one(T) - α*α + β *= one(T) - abs2(α) α = -r[k+1] for j = 1:k α -= r[k-j+1]*y[j] @@ -14,15 +14,15 @@ function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T<:BlasReal α /= β for j = 1:div(k,2) tmp = y[j] - y[j] += α*y[k-j+1] - y[k-j+1] += α*tmp + y[j] += α*conj(y[k-j+1]) + y[k-j+1] += α*conj(tmp) end - if isodd(k) y[div(k,2)+1] *= one(T) + α end + if isodd(k) y[div(k,2)+1] += α*conj(y[div(k,2)+1]) end y[k+1] = α end return y end -durbin(r::AbstractVector{T}) where {T<:BlasReal} = durbin!(r, zeros(T, length(r))) +durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r))) function levinson!(r::AbstractVector{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T<:BlasReal n = length(b) diff --git a/test/signalcorr.jl b/test/signalcorr.jl index 3234ce060..af7c33efe 100644 --- a/test/signalcorr.jl +++ b/test/signalcorr.jl @@ -114,34 +114,17 @@ c22 = crosscor(x2, x2) ## pacf -x_pacf = [-2.133252557240862 - .1775816414485478 - -.6264517920318317 - -.8809042583216906 - .09251017186697393 - -.9271887119115569 - 3.355819743178915 - -.2834039258495755 - .5354280026977677 - .39182285417742585] - -rpacfr = [-0.2181581223814188 - 0.1950153168287112 - 0.1443158046061393 - -0.1997912294497786 - 1.686452299291934 - -1.446624527965185] - -@test pacf(real.(x_pacf), 1:4) ≈ rpacfr[1:4] - -rpacfy = [-0.2211730116688735 - 0.1896833143080208 - 0.1118570207337193 - -0.1750206698354201 - 0.04490978543108597 - -0.3715806391613901 - -0.1743836518336153 - 0.07574637899230210 - 0.04970724340234131] - -@test pacf(real.(x_pacf), 1:4, method=:yulewalker) ≈ rpacfy[1:4] + +pacfr = [-1.598495044296996e-03 - 2.915104118351207e-01im + -5.560162016912027e-01 + 2.950837739894279e-01im + -2.547001916363494e-02 + 2.326084658014266e-01im + -5.427443903358727e-01 + 3.146715147305132e-01im]; + +@test pacf(x1, 1:4) ≈ pacfr[1:4] + +pacfy = [-1.598495044296996e-03 - 2.915104118351207e-01im + -5.560162016912027e-01 + 2.950837739894279e-01im + -2.547001916363494e-02 + 2.326084658014266e-01im + -5.427443903358727e-01 + 3.146715147305132e-01im]; + +@test pacf(x1, 1:4, method=:yulewalker) ≈ pacfy From 710317841fcff14f52d0b0cbf9dc4889b6254301 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Tue, 30 Nov 2021 12:21:06 +0300 Subject: [PATCH 04/15] Optimized the use of Durbin algorithm --- src/signalcorr.jl | 2 +- src/toeplitzsolvers.jl | 10 +++++---- test/signalcorr.jl | 47 ++++++++++++++++++++++++++++-------------- 3 files changed, 39 insertions(+), 20 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 51e4cb1d0..3c69f785f 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -547,7 +547,7 @@ crosscor(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) = # ####################################### -function pacf_regress!(r::AbstractMatrix X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer) +function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector, mk::Integer) lx = size(X, 1) tmpX = ones(eltype(X), lx, mk + 1) for j = 1 : size(X,2) diff --git a/src/toeplitzsolvers.jl b/src/toeplitzsolvers.jl index 28ff92423..8a8278224 100644 --- a/src/toeplitzsolvers.jl +++ b/src/toeplitzsolvers.jl @@ -1,8 +1,9 @@ # Symmetric Toeplitz solver -function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T +function durbin!(r::AbstractVector{T}, p::AbstractVector{T}, y::AbstractVector{T}) where T n = length(r) - n <= length(y) || throw(DimensionMismatch("Auxiliary vector cannot be shorter than data vector")) + n <= length(p) || n <= length(y) || throw(DimensionMismatch("Auxiliary vectors cannot be shorter than data vector")) y[1] = -r[1] + p[1] = -r[1] β = one(T) α = -r[1] for k = 1:n-1 @@ -19,10 +20,11 @@ function durbin!(r::AbstractVector{T}, y::AbstractVector{T}) where T end if isodd(k) y[div(k,2)+1] += α*conj(y[div(k,2)+1]) end y[k+1] = α + p[k+1] = α end - return y + return p end -durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r))) +durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r)), zeros(T, length(r))) function levinson!(r::AbstractVector{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T<:BlasReal n = length(b) diff --git a/test/signalcorr.jl b/test/signalcorr.jl index af7c33efe..b01220a66 100644 --- a/test/signalcorr.jl +++ b/test/signalcorr.jl @@ -113,18 +113,35 @@ c22 = crosscor(x2, x2) crosscor([34566.5345, 3466.4566], Float16[1, 10]) -## pacf - -pacfr = [-1.598495044296996e-03 - 2.915104118351207e-01im - -5.560162016912027e-01 + 2.950837739894279e-01im - -2.547001916363494e-02 + 2.326084658014266e-01im - -5.427443903358727e-01 + 3.146715147305132e-01im]; - -@test pacf(x1, 1:4) ≈ pacfr[1:4] - -pacfy = [-1.598495044296996e-03 - 2.915104118351207e-01im - -5.560162016912027e-01 + 2.950837739894279e-01im - -2.547001916363494e-02 + 2.326084658014266e-01im - -5.427443903358727e-01 + 3.146715147305132e-01im]; - -@test pacf(x1, 1:4, method=:yulewalker) ≈ pacfy +## pacf least squares +pacf_ls = [-1.598495044296996e-03 - 2.915104118351207e-01im + -5.560162016912027e-01 + 2.950837739894279e-01im + -2.547001916363494e-02 + 2.326084658014266e-01im + -5.427443903358727e-01 + 3.146715147305132e-01im] + +@test pacf(x1, 1:4) ≈ pacf_ls[1:4] + +## pacf Yule-Walker + +function yulewalker_qr(v::AbstractVector) + A = toeplitz(v) + b = v[2:end] + x =-A\b +end +function toeplitz(v::AbstractVector{T}) where T + N=length(v) + A = zeros(T, N - 1, N - 1) + for n in 1:N-1 + A[n, n + 1:end] = conj(v[2:N-n]) + A[n, 1:n] = reverse(v[1:n]) + end + return A +end +# durbin solver +acf = autocor(x1) +p = [yulewalker_qr(acf[1:n])[n-1] for n in 2:length(acf)] +@test p ≈ StatsBase.durbin(acf[2:end]) + +pacfy = []; + +@test pacf(x1, 1:4, method=:yulewalker) ≈ -p[1:4] From d32b14d3d8a6f5db76e5bc276ee3a41534af0f02 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Tue, 30 Nov 2021 13:18:19 +0300 Subject: [PATCH 05/15] Reduced allocations in autocov!() --- src/signalcorr.jl | 102 +++++++++++++++++++--------------------------- 1 file changed, 43 insertions(+), 59 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 3c69f785f..83a96fa61 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -60,33 +60,31 @@ where each column in the result will correspond to a column in `x`. The output is not normalized. See [`autocor!`](@ref) for a method with normalization. """ -function autocov!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) +function autocov!( + r::AbstractVector, + x::AbstractVector, + lags::AbstractVector{<:Integer}, + z=zeros(typeof(zero(eltype(x)) / 1), length(x)); + demean::Bool=true + ) lx = length(x) m = length(lags) length(r) == m || throw(DimensionMismatch()) check_lags(lx, lags) T = typeof(zero(eltype(x)) / 1) - z::Vector{T} = demean ? x .- mean(x) : x + demean ? z = x : z .= x .- mean(x) for k = 1 : m # foreach lag value r[k] = _autodot(z, lx, lags[k]) / lx end return r end -function autocov!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) - lx = size(x, 1) - ns = size(x, 2) - m = length(lags) - size(r) == (m, ns) || throw(DimensionMismatch()) - check_lags(lx, lags) - - T = typeof(zero(eltype(x)) / 1) - z = Vector{T}(undef, lx) - for j = 1 : ns - demean_col!(z, x, j, demean) - for k = 1 : m - r[k,j] = _autodot(z, lx, lags[k]) / lx - end +function autocov!( + r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true + ) + z = zeros(typeof(zero(eltype(x))/1), size(x, 1)) + for n in 1:size(x, 2) + autocov!(view(r, :, n), view(x, :, n), lags, z; demean) end return r end @@ -110,16 +108,16 @@ The output is not normalized. See [`autocor`](@ref) for a function with normaliz """ function autocov(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) - autocov!(out, x, lags; demean=demean) + autocov!(out, x, lags; demean) end -function autocov(x::AbstractMatrix, lags::AbstractVector; demean::Bool=true) +function autocov(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) - autocov!(out, x, lags; demean=demean) + autocov!(out, x, lags; demean) end autocov(x::AbstractVecOrMat; demean::Bool=true) = - autocov(x, default_autolags(size(x,1)); demean=demean) + autocov(x, default_autolags(size(x,1)); demean) ## autocor @@ -137,36 +135,27 @@ where each column in the result will correspond to a column in `x`. The output is normalized by the variance of `x`, i.e. so that the lag 0 autocorrelation is 1. See [`autocov!`](@ref) for the unnormalized form. """ -function autocor!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) - lx = length(x) - m = length(lags) - length(r) == m || throw(DimensionMismatch()) - check_lags(lx, lags) - - T = typeof(zero(eltype(x)) / 1) - z::Vector{T} = demean ? x .- mean(x) : x - zz = dot(z, z) - for k = 1 : m # foreach lag value - r[k] = _autodot(z, lx, lags[k]) / zz - end +function autocor!( + r::AbstractVector, + x::AbstractVector, + lags::AbstractVector{<:Integer}, + z=zeros(typeof(zero(eltype(x)) / 1), length(x)); + demean::Bool=true + ) + autocov!(view(r, 1:1), x, 0:0, z; demean) + zz = r[1] + autocov!(r, x, lags, z; demean) + r ./= zz return r end -function autocor!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) - lx = size(x, 1) - ns = size(x, 2) - m = length(lags) - size(r) == (m, ns) || throw(DimensionMismatch()) - check_lags(lx, lags) - +function autocor!( + r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true + ) T = typeof(zero(eltype(x))/1) - z = Vector{T}(undef, lx) - for j = 1 : ns - demean_col!(z, x, j, demean) - zz = dot(z, z) - for k = 1 : m - r[k,j] = _autodot(z, lx, lags[k]) / zz - end + z = Vector{T}(undef, size(x, 1)) + for n in 1:size(x, 2) + autocor!(view(r, :, n), view(x, :, n), lags, z; demean) end return r end @@ -191,16 +180,16 @@ autocorrelation is 1. See [`autocov`](@ref) for the unnormalized form. """ function autocor(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) - autocor!(out, x, lags; demean=demean) + autocor!(out, x, lags; demean) end function autocor(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) - autocor!(out, x, lags; demean=demean) + autocor!(out, x, lags; demean) end autocor(x::AbstractVecOrMat; demean::Bool=true) = - autocor(x, default_autolags(size(x,1)); demean=demean) + autocor(x, default_autolags(size(x,1)); demean) ####################################### @@ -211,13 +200,6 @@ autocor(x::AbstractVecOrMat; demean::Bool=true) = default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l) -function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int) - if l >= 0 - dot(x, 1:(lx-l), y, (1+l):lx) - else - dot(x, (1-l):lx, y, 1:(lx+l)) - end -end function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int) if l >= 0 dot(view(x, 1:(lx-l)), view(y, (1+l):lx)) @@ -565,13 +547,15 @@ function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector r end -function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer) where T<:Union{Float32, Float64} - tmp = Vector{T}(undef, mk) +function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector, mk::Integer) + p = Vector{eltype(X)}(undef, mk) + y = Vector{eltype(X)}(undef, mk) for j = 1 : size(X,2) acfs = autocor(X[:,j], 1:mk) + durbin!(acfs, p, y) for i = 1 : length(lags) l = lags[i] - r[i,j] = l == 0 ? 1 : l == 1 ? acfs[i] : -durbin!(view(acfs, 1:l), tmp)[l] + r[i,j] = l == 0 ? 1 : -p[l] end end end From 97d87d310ae27cf7fcbde7220361d1168eb4a34c Mon Sep 17 00:00:00 2001 From: Mikhail Kagalenko <16374215+kagalenko-m-b@users.noreply.github.com> Date: Tue, 30 Nov 2021 13:34:53 +0300 Subject: [PATCH 06/15] Update test/signalcorr.jl Co-authored-by: Moritz Schauer --- test/signalcorr.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/signalcorr.jl b/test/signalcorr.jl index b01220a66..ae8ba4cb8 100644 --- a/test/signalcorr.jl +++ b/test/signalcorr.jl @@ -126,7 +126,7 @@ pacf_ls = [-1.598495044296996e-03 - 2.915104118351207e-01im function yulewalker_qr(v::AbstractVector) A = toeplitz(v) b = v[2:end] - x =-A\b + x = -A\b end function toeplitz(v::AbstractVector{T}) where T N=length(v) From 1fff4e7bd324760be5e66bccddaea84c61db16c5 Mon Sep 17 00:00:00 2001 From: Mikhail Kagalenko <16374215+kagalenko-m-b@users.noreply.github.com> Date: Tue, 30 Nov 2021 17:20:45 +0300 Subject: [PATCH 07/15] Apply suggestions from code review Implement assorted suggestions from code review Co-authored-by: David Widmann Co-authored-by: Moritz Schauer --- src/signalcorr.jl | 14 +++++++------- src/toeplitzsolvers.jl | 4 +++- test/signalcorr.jl | 2 +- 3 files changed, 11 insertions(+), 9 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 83a96fa61..8f99b7013 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -108,16 +108,16 @@ The output is not normalized. See [`autocor`](@ref) for a function with normaliz """ function autocov(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) - autocov!(out, x, lags; demean) + autocov!(out, x, lags; demean=demean) end function autocov(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) - autocov!(out, x, lags; demean) + autocov!(out, x, lags; demean=demean) end autocov(x::AbstractVecOrMat; demean::Bool=true) = - autocov(x, default_autolags(size(x,1)); demean) + autocov(x, default_autolags(size(x,1)); demean=demean) ## autocor @@ -155,7 +155,7 @@ function autocor!( T = typeof(zero(eltype(x))/1) z = Vector{T}(undef, size(x, 1)) for n in 1:size(x, 2) - autocor!(view(r, :, n), view(x, :, n), lags, z; demean) + autocor!(view(r, :, n), view(x, :, n), lags, z; demean=demean) end return r end @@ -180,16 +180,16 @@ autocorrelation is 1. See [`autocov`](@ref) for the unnormalized form. """ function autocor(x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Vector{float(eltype(x))}(undef, length(lags)) - autocor!(out, x, lags; demean) + autocor!(out, x, lags; demean=demean) end function autocor(x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) out = Matrix{float(eltype(x))}(undef, length(lags), size(x,2)) - autocor!(out, x, lags; demean) + autocor!(out, x, lags; demean=demean) end autocor(x::AbstractVecOrMat; demean::Bool=true) = - autocor(x, default_autolags(size(x,1)); demean) + autocor(x, default_autolags(size(x,1)); demean=demean) ####################################### diff --git a/src/toeplitzsolvers.jl b/src/toeplitzsolvers.jl index 8a8278224..81be20c4c 100644 --- a/src/toeplitzsolvers.jl +++ b/src/toeplitzsolvers.jl @@ -18,7 +18,9 @@ function durbin!(r::AbstractVector{T}, p::AbstractVector{T}, y::AbstractVector{T y[j] += α*conj(y[k-j+1]) y[k-j+1] += α*conj(tmp) end - if isodd(k) y[div(k,2)+1] += α*conj(y[div(k,2)+1]) end + if isodd(k) + y[div(k,2)+1] += α*conj(y[div(k,2)+1]) + end y[k+1] = α p[k+1] = α end diff --git a/test/signalcorr.jl b/test/signalcorr.jl index ae8ba4cb8..a2b4ddcd8 100644 --- a/test/signalcorr.jl +++ b/test/signalcorr.jl @@ -132,7 +132,7 @@ function toeplitz(v::AbstractVector{T}) where T N=length(v) A = zeros(T, N - 1, N - 1) for n in 1:N-1 - A[n, n + 1:end] = conj(v[2:N-n]) + A[n, n+1:end] = conj(v[2:N-n]) A[n, 1:n] = reverse(v[1:n]) end return A From 57b0fab7971c67aa0edc18f52bdc7fa2c5874052 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Tue, 30 Nov 2021 19:57:40 +0300 Subject: [PATCH 08/15] Add documentation to durbin!() and restore tests --- src/signalcorr.jl | 14 ++++---- src/toeplitzsolvers.jl | 28 +++++++++++++-- test/signalcorr.jl | 78 +++++++++++++++++++++++++----------------- 3 files changed, 79 insertions(+), 41 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 8f99b7013..63a0739d7 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -45,7 +45,6 @@ default_autolags(lx::Int) = 0 : default_laglen(lx) _autodot(x::AbstractVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx)) - ## autocov """ autocov!(r, x, lags; demean=true) @@ -64,15 +63,15 @@ function autocov!( r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}, - z=zeros(typeof(zero(eltype(x)) / 1), length(x)); + z=Vector{typeof(zero(eltype(x)) / 1)}(undef, size(x, 1)); demean::Bool=true ) lx = length(x) m = length(lags) length(r) == m || throw(DimensionMismatch()) check_lags(lx, lags) - T = typeof(zero(eltype(x)) / 1) - demean ? z = x : z .= x .- mean(x) + z .= x + demean && z .= z .- mean(z) for k = 1 : m # foreach lag value r[k] = _autodot(z, lx, lags[k]) / lx end @@ -82,7 +81,8 @@ end function autocov!( r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true ) - z = zeros(typeof(zero(eltype(x))/1), size(x, 1)) + T = typeof(zero(eltype(x)) / 1) + z = Vector{T}(undef, size(x, 1)) for n in 1:size(x, 2) autocov!(view(r, :, n), view(x, :, n), lags, z; demean) end @@ -145,7 +145,7 @@ function autocor!( autocov!(view(r, 1:1), x, 0:0, z; demean) zz = r[1] autocov!(r, x, lags, z; demean) - r ./= zz + ldiv!(zz, r) return r end @@ -552,7 +552,7 @@ function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVec y = Vector{eltype(X)}(undef, mk) for j = 1 : size(X,2) acfs = autocor(X[:,j], 1:mk) - durbin!(acfs, p, y) + durbin!(acfs, y, p) for i = 1 : length(lags) l = lags[i] r[i,j] = l == 0 ? 1 : -p[l] diff --git a/src/toeplitzsolvers.jl b/src/toeplitzsolvers.jl index 81be20c4c..bc254c3cd 100644 --- a/src/toeplitzsolvers.jl +++ b/src/toeplitzsolvers.jl @@ -1,5 +1,27 @@ -# Symmetric Toeplitz solver -function durbin!(r::AbstractVector{T}, p::AbstractVector{T}, y::AbstractVector{T}) where T +""" + durbin!(r::AbstractVector, y::AbstractVector, p::AbstractVector) -> Vector + +Solve Yule-Walker equations using Durbin algorithm. + +Solution of NxN system is obtained iteratively, by solving 1x1, 2x2, ... +in succession. For use in computing partial autocorrelation matrix, +return the last elements of the successive solutions in vector p[]. + +The section 4.7 of Golub,VanLoan "Matrix Computations," 4th ed., +discusses the algorithm in detail. + +# Arguments +- `r::AbstractVector`: first column of a herimitian positive definite + Toeplitz matrix, excluding the diagonal element (equal to one). +- `y::AbstractVector`: work vector for solution, should have the same length + as r[] +- `p::AbstractVector`: last elements of the successive solutions, should + have the same length as r[] + +# Returns +- `y::AbstractVector`: same as the second argument +""" +function durbin!(r::AbstractVector{T}, y::AbstractVector{T}, p::AbstractVector{T}) where T n = length(r) n <= length(p) || n <= length(y) || throw(DimensionMismatch("Auxiliary vectors cannot be shorter than data vector")) y[1] = -r[1] @@ -24,7 +46,7 @@ function durbin!(r::AbstractVector{T}, p::AbstractVector{T}, y::AbstractVector{T y[k+1] = α p[k+1] = α end - return p + return y end durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r)), zeros(T, length(r))) diff --git a/test/signalcorr.jl b/test/signalcorr.jl index a2b4ddcd8..566a8280b 100644 --- a/test/signalcorr.jl +++ b/test/signalcorr.jl @@ -22,39 +22,46 @@ x = [0.5316732188954801 + 0.015032026010878084im -0.8051436916816284 - 0.7302380 x1 = view(x, :, 1) x2 = view(x, :, 2) - +cmplx_x = convert(AbstractMatrix{Complex}, x) +cmplx_x1 = convert(AbstractVector{Complex}, x1) +cmplx_x2 = convert(AbstractVector{Complex}, x2) # autocov & autocorr @test autocov([1:5;]) ≈ [2.0, 0.8, -0.2, -0.8, -0.8] @test autocor([1, 2, 3, 4, 5]) ≈ [1.0, 0.4, -0.1, -0.4, -0.4] -racovx1 = [0.755284179631112 + 0.0im, - -0.005333112170365584 - 0.18633291805458002im, - -0.28638133506011604 + 0.1397225175604478im, - -0.05476759020476188 + 0.12991227531087618im, - -0.00499902418309206 + 0.023463421186162473im, - 0.11854989361702409 - 0.0028772373657594066im, - -0.005123812003979093 - 0.08041878599400383im, - -0.14090692583679154 + 0.035230042290069444im, - 0.039899180711674635 + 0.004918236900383659im, - -0.03857936468514856 - 0.04063999068714769im] - -@test autocov(x1) ≈ racovx1 +acovx1 = [0.755284179631112 + 0.0im, + -0.005333112170365584 - 0.18633291805458002im, + -0.28638133506011604 + 0.1397225175604478im, + -0.05476759020476188 + 0.12991227531087618im, + -0.00499902418309206 + 0.023463421186162473im, + 0.11854989361702409 - 0.0028772373657594066im, + -0.005123812003979093 - 0.08041878599400383im, + -0.14090692583679154 + 0.035230042290069444im, + 0.039899180711674635 + 0.004918236900383659im, + -0.03857936468514856 - 0.04063999068714769im] + +@test autocov(x1) ≈ acovx1 +@test autocov(cmplx_x1) ≈ acovx1 @test autocov(x) ≈ [autocov(x1) autocov(x2)] - -racorx1 = [1.0 + 0.0im, - -0.007061066965510023 - 0.24670570770539213im, - -0.3791703080554228 + 0.18499330626611243im, - -0.07251256107537019 + 0.17200449686941222im, - -0.006618732813301651 + 0.031065686027770673im, - 0.1569606471499575 - 0.0038094765432061303im, - -0.00678395250709688 - 0.106474871528861im, - -0.18656146869859214 + 0.04664475073114351im, - 0.05282671316002114 + 0.0065117700502952056im, - -0.051079270194684986 - 0.0538075492419246im] - -@test autocor(x1) ≈ racorx1 +@test autocov(cmplx_x) ≈ [autocov(cmplx_x1) autocov(cmplx_x2)] + +acorx1 = [1.0 + 0.0im, + -0.007061066965510023 - 0.24670570770539213im, + -0.3791703080554228 + 0.18499330626611243im, + -0.07251256107537019 + 0.17200449686941222im, + -0.006618732813301651 + 0.031065686027770673im, + 0.1569606471499575 - 0.0038094765432061303im, + -0.00678395250709688 - 0.106474871528861im, + -0.18656146869859214 + 0.04664475073114351im, + 0.05282671316002114 + 0.0065117700502952056im, + -0.051079270194684986 - 0.0538075492419246im] + +@test autocor(x1) ≈ acorx1 +@test autocor(cmplx_x1) ≈ acorx1 @test autocor(x) ≈ [autocor(x1) autocor(x2)] +@test autocor(cmplx_x) ≈ [autocor(cmplx_x1) autocor(cmplx_x2)] + # crosscov & crosscor @@ -75,10 +82,14 @@ c11 = crosscov(x1, x1) c12 = crosscov(x1, x2) c21 = crosscov(x2, x1) c22 = crosscov(x2, x2) +@test crosscov(cmplx_x1, cmplx_x2) ≈ c12 @test crosscov(x, x1) ≈ [c11 c21] +@test crosscov(cmplx_x, cmplx_x1) ≈ [c11 c21] @test crosscov(x1, x) ≈ [c11 c12] +@test crosscov(cmplx_x1, cmplx_x) ≈ [c11 c12] @test crosscov(x, x) ≈ cat([c11 c21], [c12 c22], dims=3) +@test crosscov(cmplx_x, cmplx_x) ≈ cat([c11 c21], [c12 c22], dims=3) # issue #805: avoid converting one input to the other's eltype @test crosscov([34566.5345, 3466.4566], Float16[1, 10]) ≈ @@ -102,10 +113,14 @@ c11 = crosscor(x1, x1) c12 = crosscor(x1, x2) c21 = crosscor(x2, x1) c22 = crosscor(x2, x2) +@test crosscor(cmplx_x1, cmplx_x2) ≈ c12 @test crosscor(x, x1) ≈ [c11 c21] +@test crosscor(cmplx_x, cmplx_x1) ≈ [c11 c21] @test crosscor(x1, x) ≈ [c11 c12] +@test crosscor(cmplx_x1, cmplx_x) ≈ [c11 c12] @test crosscor(x, x) ≈ cat([c11 c21], [c12 c22], dims=3) +@test crosscor(cmplx_x, cmplx_x) ≈ cat([c11 c21], [c12 c22], dims=3) # issue #805: avoid converting one input to the other's eltype @test crosscor([34566.5345, 3466.4566], Float16[1, 10]) ≈ @@ -139,9 +154,10 @@ function toeplitz(v::AbstractVector{T}) where T end # durbin solver acf = autocor(x1) -p = [yulewalker_qr(acf[1:n])[n-1] for n in 2:length(acf)] -@test p ≈ StatsBase.durbin(acf[2:end]) - -pacfy = []; +p_qr = [yulewalker_qr(acf[1:n])[n-1] for n in 2:length(acf)] +y = similar(acf[2:end]) +pd = similar(acf[2:end]) +StatsBase.durbin!(acf[2:end], y, pd) +@test p_qr ≈ pd -@test pacf(x1, 1:4, method=:yulewalker) ≈ -p[1:4] +@test pacf(x1, 1:4, method=:yulewalker) ≈ -p_qr[1:4] From c06c4d7fede00b763c4ac18c5ec2882e0e22db68 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Wed, 1 Dec 2021 13:11:57 +0300 Subject: [PATCH 09/15] Simplify subtracting mean, as suggested by reviewer --- src/signalcorr.jl | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 63a0739d7..33b401e73 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -70,8 +70,7 @@ function autocov!( m = length(lags) length(r) == m || throw(DimensionMismatch()) check_lags(lx, lags) - z .= x - demean && z .= z .- mean(z) + demean ? z .= x .- mean(x) : copyto!(z, x) for k = 1 : m # foreach lag value r[k] = _autodot(z, lx, lags[k]) / lx end @@ -529,7 +528,7 @@ crosscor(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) = # ####################################### -function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector, mk::Integer) +function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer) lx = size(X, 1) tmpX = ones(eltype(X), lx, mk + 1) for j = 1 : size(X,2) @@ -547,7 +546,7 @@ function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector r end -function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector, mk::Integer) +function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer) p = Vector{eltype(X)}(undef, mk) y = Vector{eltype(X)}(undef, mk) for j = 1 : size(X,2) @@ -572,7 +571,7 @@ using the Yule-Walker equations. `r` must be a matrix of size `(length(lags), size(x, 2))`. """ -function pacf!(r::AbstractMatrix, X::AbstractMatrix, lags::IntegerVector; method::Symbol=:regression) +function pacf!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}; method::Symbol=:regression) lx = size(X, 1) m = length(lags) minlag, maxlag = extrema(lags) @@ -603,11 +602,11 @@ If `x` is a vector, return a vector of the same length as `lags`. If `x` is a matrix, return a matrix of size `(length(lags), size(x, 2))`, where each column in the result corresponds to a column in `x`. """ -function pacf(X::AbstractMatrix, lags::IntegerVector; method::Symbol=:regression) +function pacf(X::AbstractMatrix, lags::AbstractVector{<:Integer}; method::Symbol=:regression) out = Matrix{float(eltype(X))}(undef, length(lags), size(X,2)) pacf!(out, float(X), lags; method=method) end -function pacf(x::AbstractVector, lags::IntegerVector; method::Symbol=:regression) +function pacf(x::AbstractVector, lags::AbstractVector{<:Integer}; method::Symbol=:regression) vec(pacf(reshape(x, length(x), 1), lags, method=method)) end From dac79be8743371a038f154bad5e46f9c264096bf Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Thu, 2 Dec 2021 09:54:28 +0300 Subject: [PATCH 10/15] Moved new tests with complex data into a separate file --- test/{signalcorr.jl => signalcorr_complex.jl} | 0 test/signalcorr_real.jl | 136 ++++++++++++++++++ 2 files changed, 136 insertions(+) rename test/{signalcorr.jl => signalcorr_complex.jl} (100%) create mode 100644 test/signalcorr_real.jl diff --git a/test/signalcorr.jl b/test/signalcorr_complex.jl similarity index 100% rename from test/signalcorr.jl rename to test/signalcorr_complex.jl diff --git a/test/signalcorr_real.jl b/test/signalcorr_real.jl new file mode 100644 index 000000000..bce1c83a3 --- /dev/null +++ b/test/signalcorr_real.jl @@ -0,0 +1,136 @@ +# Test corr.jl +# +# Many test cases are migrated from test/01.jl in the old version +# The reference results are generated from R. +# + +using StatsBase +using Test + +# random data for testing + +x = [-2.133252557240862 -.7445937365828654; + .1775816414485478 -.5834801838041446; + -.6264517920318317 -.68444205333293; + -.8809042583216906 .9071671734302398; + .09251017186697393 -1.0404476733379926; + -.9271887119115569 -.620728578941385; + 3.355819743178915 -.8325051361909978; + -.2834039258495755 -.22394811874731657; + .5354280026977677 .7481337671592626; + .39182285417742585 .3085762550821047] + +x1 = view(x, :, 1) +x2 = view(x, :, 2) +realx = convert(AbstractMatrix{Real}, x) +realx1 = convert(AbstractVector{Real}, x1) +realx2 = convert(AbstractVector{Real}, x2) + +# autocov & autocorr + +@test autocov([1:5;]) ≈ [2.0, 0.8, -0.2, -0.8, -0.8] +@test autocor([1, 2, 3, 4, 5]) ≈ [1.0, 0.4, -0.1, -0.4, -0.4] + +racovx1 = [1.839214242630635709475, + -0.406784553146903871124, + 0.421772254824993531042, + 0.035874943792884653182, + -0.255679775928512320604, + 0.231154400105831353551, + -0.787016960267425180753, + 0.039909287349160660341, + -0.110149697877911914579, + -0.088687020167434751916] + +@test autocov(x1) ≈ racovx1 +@test autocov(realx1) ≈ racovx1 +@test autocov(x) ≈ [autocov(x1) autocov(x2)] +@test autocov(realx) ≈ [autocov(realx1) autocov(realx2)] + +racorx1 = [0.999999999999999888978, + -0.221173011668873431557, + 0.229321981664153962122, + 0.019505581764945757045, + -0.139015765538446717242, + 0.125681062460244019618, + -0.427909344123907742219, + 0.021699096507690283225, + -0.059889541590524189574, + -0.048220059475281865091] + +@test autocor(x1) ≈ racorx1 +@test autocor(realx1) ≈ racorx1 +@test autocor(x) ≈ [autocor(x1) autocor(x2)] +@test autocor(realx) ≈ [autocor(realx1) autocor(realx2)] + + +# crosscov & crosscor + +rcov0 = [0.320000000000000006661, + -0.319999999999999951150, + 0.080000000000000029421, + -0.479999999999999982236, + 0.000000000000000000000, + 0.479999999999999982236, + -0.080000000000000029421, + 0.319999999999999951150, + -0.320000000000000006661] + +@test crosscov([1, 2, 3, 4, 5], [1, -1, 1, -1, 1]) ≈ rcov0 +@test crosscov([1:5;], [1:5;]) ≈ [-0.8, -0.8, -0.2, 0.8, 2.0, 0.8, -0.2, -0.8, -0.8] + +c11 = crosscov(x1, x1) +c12 = crosscov(x1, x2) +c21 = crosscov(x2, x1) +c22 = crosscov(x2, x2) +@test crosscov(realx1, realx2) ≈ c12 + +@test crosscov(x, x1) ≈ [c11 c21] +@test crosscov(realx, realx1) ≈ [c11 c21] +@test crosscov(x1, x) ≈ [c11 c12] +@test crosscov(realx1, realx) ≈ [c11 c12] +@test crosscov(x, x) ≈ cat([c11 c21], [c12 c22], dims=3) +@test crosscov(realx, realx) ≈ cat([c11 c21], [c12 c22], dims=3) + +rcor0 = [0.230940107675850, + -0.230940107675850, + 0.057735026918963, + -0.346410161513775, + 0.000000000000000, + 0.346410161513775, + -0.057735026918963, + 0.230940107675850, + -0.230940107675850] + +@test crosscor([1, 2, 3, 4, 5], [1, -1, 1, -1, 1]) ≈ rcor0 +@test crosscor([1:5;], [1:5;]) ≈ [-0.4, -0.4, -0.1, 0.4, 1.0, 0.4, -0.1, -0.4, -0.4] + +c11 = crosscor(x1, x1) +c12 = crosscor(x1, x2) +c21 = crosscor(x2, x1) +c22 = crosscor(x2, x2) +@test crosscor(realx1, realx2) ≈ c12 + +@test crosscor(x, x1) ≈ [c11 c21] +@test crosscor(realx, realx1) ≈ [c11 c21] +@test crosscor(x1, x) ≈ [c11 c12] +@test crosscor(realx1, realx) ≈ [c11 c12] +@test crosscor(x, x) ≈ cat([c11 c21], [c12 c22], dims=3) +@test crosscor(realx, realx) ≈ cat([c11 c21], [c12 c22], dims=3) + + +## pacf + +rpacfr = [-0.218158122381419, + 0.195015316828711, + 0.144315804606139, + -0.199791229449779] + +@test pacf(x[:,1], 1:4) ≈ rpacfr + +rpacfy = [-0.221173011668873, + 0.189683314308021, + 0.111857020733719, + -0.175020669835420] + +@test pacf(x[:,1], 1:4, method=:yulewalker) ≈ rpacfy From 3384cb153e9f483ad3913fd20e0b0cceaf89fc8a Mon Sep 17 00:00:00 2001 From: Mikhail Kagalenko <16374215+kagalenko-m-b@users.noreply.github.com> Date: Fri, 3 Dec 2021 16:49:06 +0300 Subject: [PATCH 11/15] Apply suggestions from code review Co-authored-by: David Widmann --- src/signalcorr.jl | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 33b401e73..4cbe98951 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -63,23 +63,23 @@ function autocov!( r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}, - z=Vector{typeof(zero(eltype(x)) / 1)}(undef, size(x, 1)); + z::Vector=Vector{typeof(zero(eltype(x)) / 1)}(undef, length(x)); demean::Bool=true - ) +) lx = length(x) m = length(lags) - length(r) == m || throw(DimensionMismatch()) + length(r) == length(z) == m || throw(DimensionMismatch()) check_lags(lx, lags) demean ? z .= x .- mean(x) : copyto!(z, x) - for k = 1 : m # foreach lag value - r[k] = _autodot(z, lx, lags[k]) / lx + @inbounds for (k, lags_k) in zip(eachindex(r), lags) # foreach lag value + r[k] = _autodot(z, lx, lags_k) / lx end return r end function autocov!( r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true - ) +) T = typeof(zero(eltype(x)) / 1) z = Vector{T}(undef, size(x, 1)) for n in 1:size(x, 2) From 9f2c943f0a95cc86823cedb240935607c8c6f192 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Sat, 4 Dec 2021 11:28:05 +0300 Subject: [PATCH 12/15] Limit changes to more general argument types --- src/signalcorr.jl | 91 +++++++++++++++++++++++++++++------------------ 1 file changed, 56 insertions(+), 35 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 4cbe98951..955f83884 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -43,8 +43,10 @@ end default_autolags(lx::Int) = 0 : default_laglen(lx) +_autodot(x::AbstractVector{<:RealFP}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx) _autodot(x::AbstractVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx)) + ## autocov """ autocov!(r, x, lags; demean=true) @@ -59,31 +61,34 @@ where each column in the result will correspond to a column in `x`. The output is not normalized. See [`autocor!`](@ref) for a method with normalization. """ -function autocov!( - r::AbstractVector, - x::AbstractVector, - lags::AbstractVector{<:Integer}, - z::Vector=Vector{typeof(zero(eltype(x)) / 1)}(undef, length(x)); - demean::Bool=true -) +function autocov!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = length(x) m = length(lags) - length(r) == length(z) == m || throw(DimensionMismatch()) + length(r) == m || throw(DimensionMismatch()) check_lags(lx, lags) - demean ? z .= x .- mean(x) : copyto!(z, x) - @inbounds for (k, lags_k) in zip(eachindex(r), lags) # foreach lag value - r[k] = _autodot(z, lx, lags_k) / lx + + T = typeof(zero(eltype(x)) / 1) + z::Vector{T} = demean ? x .- mean(x) : x + for k = 1 : m # foreach lag value + r[k] = _autodot(z, lx, lags[k]) / lx end return r end -function autocov!( - r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true -) +function autocov!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) + lx = size(x, 1) + ns = size(x, 2) + m = length(lags) + size(r) == (m, ns) || throw(DimensionMismatch()) + check_lags(lx, lags) + T = typeof(zero(eltype(x)) / 1) - z = Vector{T}(undef, size(x, 1)) - for n in 1:size(x, 2) - autocov!(view(r, :, n), view(x, :, n), lags, z; demean) + z = Vector{T}(undef, lx) + for j = 1 : ns + demean_col!(z, x, j, demean) + for k = 1 : m + r[k,j] = _autodot(z, lx, lags[k]) / lx + end end return r end @@ -134,27 +139,36 @@ where each column in the result will correspond to a column in `x`. The output is normalized by the variance of `x`, i.e. so that the lag 0 autocorrelation is 1. See [`autocov!`](@ref) for the unnormalized form. """ -function autocor!( - r::AbstractVector, - x::AbstractVector, - lags::AbstractVector{<:Integer}, - z=zeros(typeof(zero(eltype(x)) / 1), length(x)); - demean::Bool=true - ) - autocov!(view(r, 1:1), x, 0:0, z; demean) - zz = r[1] - autocov!(r, x, lags, z; demean) - ldiv!(zz, r) +function autocor!(r::AbstractVector, x::AbstractVector, lags::AbstractVector{<:Integer}; demean::Bool=true) + lx = length(x) + m = length(lags) + length(r) == m || throw(DimensionMismatch()) + check_lags(lx, lags) + + T = typeof(zero(eltype(x)) / 1) + z::Vector{T} = demean ? x .- mean(x) : x + zz = dot(z, z) + for k = 1 : m # foreach lag value + r[k] = _autodot(z, lx, lags[k]) / zz + end return r end -function autocor!( - r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true - ) - T = typeof(zero(eltype(x))/1) - z = Vector{T}(undef, size(x, 1)) - for n in 1:size(x, 2) - autocor!(view(r, :, n), view(x, :, n), lags, z; demean=demean) +function autocor!(r::AbstractMatrix, x::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) + lx = size(x, 1) + ns = size(x, 2) + m = length(lags) + size(r) == (m, ns) || throw(DimensionMismatch()) + check_lags(lx, lags) + + T = typeof(zero(eltype(x)) / 1) + z = Vector{T}(undef, lx) + for j = 1 : ns + demean_col!(z, x, j, demean) + zz = dot(z, z) + for k = 1 : m + r[k,j] = _autodot(z, lx, lags[k]) / zz + end end return r end @@ -199,6 +213,13 @@ autocor(x::AbstractVecOrMat; demean::Bool=true) = default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l) +function _crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:RealFP} + if l >= 0 + dot(x, 1:(lx-l), y, (1+l):lx) + else + dot(x, (1-l):lx, y, 1:(lx+l)) + end +end function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int) if l >= 0 dot(view(x, 1:(lx-l)), view(y, (1+l):lx)) From cb6ea6491944e87cde09ab1ea9a085ed5a70136f Mon Sep 17 00:00:00 2001 From: kagalenko-m-b Date: Sun, 5 Dec 2021 18:48:36 +0300 Subject: [PATCH 13/15] Apply suggestions from code review --- src/signalcorr.jl | 4 ++-- src/toeplitzsolvers.jl | 19 +++++++++---------- test/runtests.jl | 1 + test/{signalcorr_complex.jl => signalcorr.jl} | 0 4 files changed, 12 insertions(+), 12 deletions(-) rename test/{signalcorr_complex.jl => signalcorr.jl} (100%) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 955f83884..8a289538c 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -302,7 +302,7 @@ function crosscov!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags return r end -function crosscov!(r::AbstractArray{U,3} where U, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov!(r::AbstractArray{<:Any}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) @@ -458,7 +458,7 @@ function crosscor!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags return r end -function crosscor!(r::AbstractArray{U,3} where U, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor!(r::AbstractArray{<:Any}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) diff --git a/src/toeplitzsolvers.jl b/src/toeplitzsolvers.jl index bc254c3cd..16a905a6c 100644 --- a/src/toeplitzsolvers.jl +++ b/src/toeplitzsolvers.jl @@ -3,25 +3,24 @@ Solve Yule-Walker equations using Durbin algorithm. -Solution of NxN system is obtained iteratively, by solving 1x1, 2x2, ... +Solution of N×N system is obtained iteratively, by solving 1×1, 2×2, ... in succession. For use in computing partial autocorrelation matrix, -return the last elements of the successive solutions in vector p[]. +return the last elements of the successive solutions in vector `p`. -The section 4.7 of Golub,VanLoan "Matrix Computations," 4th ed., -discusses the algorithm in detail. +Reference: Golub, G. H., and C. F. Van Loan. "Matrix computations 4th edition the johns hopkins university press." Baltimore, MD (2013), section 4.7 # Arguments -- `r::AbstractVector`: first column of a herimitian positive definite +- `r::AbstractVector`: first column of a Herimitian positive definite Toeplitz matrix, excluding the diagonal element (equal to one). - `y::AbstractVector`: work vector for solution, should have the same length - as r[] + as `r` - `p::AbstractVector`: last elements of the successive solutions, should - have the same length as r[] + have the same length as `r` # Returns -- `y::AbstractVector`: same as the second argument +- `y::AbstractVector`: return the solution vector (same as the second argument) """ -function durbin!(r::AbstractVector{T}, y::AbstractVector{T}, p::AbstractVector{T}) where T +function durbin!(r::AbstractVector{T}, y::AbstractVector{T}, p::AbstractVector{T}) where {T} n = length(r) n <= length(p) || n <= length(y) || throw(DimensionMismatch("Auxiliary vectors cannot be shorter than data vector")) y[1] = -r[1] @@ -48,7 +47,7 @@ function durbin!(r::AbstractVector{T}, y::AbstractVector{T}, p::AbstractVector{T end return y end -durbin(r::AbstractVector{T}) where T = durbin!(r, zeros(T, length(r)), zeros(T, length(r))) +durbin(r::AbstractVector{T}) where {T} = durbin!(r, zeros(T, length(r)), zeros(T, length(r))) function levinson!(r::AbstractVector{T}, b::AbstractVector{T}, x::AbstractVector{T}) where T<:BlasReal n = length(b) diff --git a/test/runtests.jl b/test/runtests.jl index 94986d240..1caa2269a 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,6 +15,7 @@ tests = ["ambiguous", "hist", "rankcorr", "signalcorr", + "signalcorr_real", "misc", "pairwise", "reliability", diff --git a/test/signalcorr_complex.jl b/test/signalcorr.jl similarity index 100% rename from test/signalcorr_complex.jl rename to test/signalcorr.jl From c135edf096023faa7e01621714e4693db558f38a Mon Sep 17 00:00:00 2001 From: kagalenko-m-b <16374215+kagalenko-m-b@users.noreply.github.com> Date: Thu, 16 Dec 2021 11:58:56 +0300 Subject: [PATCH 14/15] Correct the documentation to reflect more general input type --- src/signalcorr.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/signalcorr.jl b/src/signalcorr.jl index 8a289538c..c419137e6 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -233,7 +233,7 @@ end """ crosscov!(r, x, y, lags; demean=true) -Compute the cross covariance function (CCF) between real-valued vectors or matrices +Compute the cross covariance function (CCF) between vectors or matrices `x` and `y` at `lags` and store the result in `r`. `demean` specifies whether the respective means of `x` and `y` should be subtracted from them before computing their CCF. @@ -343,7 +343,7 @@ end """ crosscov(x, y, [lags]; demean=true) -Compute the cross covariance function (CCF) between real-valued vectors or +Compute the cross covariance function (CCF) between vectors or matrices `x` and `y`, optionally specifying the `lags`. `demean` specifies whether the respective means of `x` and `y` should be subtracted from them before computing their CCF. @@ -384,7 +384,7 @@ crosscov(x::AbstractVecOrMat, y::AbstractVecOrMat; demean::Bool=true) = """ crosscor!(r, x, y, lags; demean=true) -Compute the cross correlation between real-valued vectors or matrices `x` and `y` at +Compute the cross correlation between vectors or matrices `x` and `y` at `lags` and store the result in `r`. `demean` specifies whether the respective means of `x` and `y` should be subtracted from them before computing their cross correlation. @@ -504,7 +504,7 @@ end """ crosscor(x, y, [lags]; demean=true) -Compute the cross correlation between real-valued vectors or matrices `x` and `y`, +Compute the cross correlation between vectors or matrices `x` and `y`, optionally specifying the `lags`. `demean` specifies whether the respective means of `x` and `y` should be subtracted from them before computing their cross correlation. @@ -613,7 +613,7 @@ end """ pacf(X, lags; method=:regression) -Compute the partial autocorrelation function (PACF) of a real-valued vector +Compute the partial autocorrelation function (PACF) of a vector or matrix `X` at `lags`. `method` designates the estimation method. Recognized values are `:regression`, which computes the partial autocorrelations via successive regression models, and `:yulewalker`, which computes the partial autocorrelations From 4f3021157ad8ff2339d55393e210ffd474d57d34 Mon Sep 17 00:00:00 2001 From: kagalenko-m-b Date: Sun, 29 May 2022 17:32:40 +0300 Subject: [PATCH 15/15] Remove unneeded stuff from tests --- src/signalcorr.jl | 13 +++++++------ test/rcall_test.jl | 29 ----------------------------- test/signalcorr.jl | 4 +++- 3 files changed, 10 insertions(+), 36 deletions(-) delete mode 100644 test/rcall_test.jl diff --git a/src/signalcorr.jl b/src/signalcorr.jl index c419137e6..99dd0385f 100644 --- a/src/signalcorr.jl +++ b/src/signalcorr.jl @@ -43,7 +43,7 @@ end default_autolags(lx::Int) = 0 : default_laglen(lx) -_autodot(x::AbstractVector{<:RealFP}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx) +_autodot(x::AbstractVector{<:Union{Float32, Float64}}, lx::Int, l::Int) = dot(x, 1:(lx-l), x, (1+l):lx) _autodot(x::AbstractVector, lx::Int, l::Int) = dot(view(x, 1:(lx-l)), view(x, (1+l):lx)) @@ -213,7 +213,7 @@ autocor(x::AbstractVecOrMat; demean::Bool=true) = default_crosslags(lx::Int) = (l=default_laglen(lx); -l:l) -function _crossdot(x::AbstractVector{T}, y::AbstractVector{T}, lx::Int, l::Int) where {T<:RealFP} +function _crossdot(x::AbstractVector, y::AbstractVector, lx::Int, l::Int) if l >= 0 dot(x, 1:(lx-l), y, (1+l):lx) else @@ -302,7 +302,7 @@ function crosscov!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags return r end -function crosscov!(r::AbstractArray{<:Any}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscov!(r::AbstractArray{<:Any,3}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) @@ -458,7 +458,7 @@ function crosscor!(r::AbstractMatrix, x::AbstractVector, y::AbstractMatrix, lags return r end -function crosscor!(r::AbstractArray{<:Any}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) +function crosscor!(r::AbstractArray{<:Any, 3}, x::AbstractMatrix, y::AbstractMatrix, lags::AbstractVector{<:Integer}; demean::Bool=true) lx = size(x, 1) nx = size(x, 2) ny = size(y, 2) @@ -568,8 +568,9 @@ function pacf_regress!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVecto end function pacf_yulewalker!(r::AbstractMatrix, X::AbstractMatrix, lags::AbstractVector{<:Integer}, mk::Integer) - p = Vector{eltype(X)}(undef, mk) - y = Vector{eltype(X)}(undef, mk) + T = eltype(X) + p = Vector{T}(undef, mk) + y = Vector{T}(undef, mk) for j = 1 : size(X,2) acfs = autocor(X[:,j], 1:mk) durbin!(acfs, y, p) diff --git a/test/rcall_test.jl b/test/rcall_test.jl deleted file mode 100644 index efa0c1bfb..000000000 --- a/test/rcall_test.jl +++ /dev/null @@ -1,29 +0,0 @@ -using RCall - -_r_ccf(x, y, lm) = rcopy(rcall(:ccf, x, y, plot=false, lag=lm, type="covariance"))[:acf] -_r_acf(x, lm) = rcopy(rcall(:acf, x, plot=false, lag=lm, type="covariance"))[:acf] -_r_pacf(x, lm) = rcopy(rcall(:acf, x, plot=false, lag=lm))[:acf] - -function _ccf(is_auto::Bool, x::AbstractVector{<:Real}...) - lx = size(x[1], 1) - lag_max = min(lx - 1, round(Int, 10*log10(lx))) - cor = reverse(dropdims(_r_ccf(x[1], x[end], lag_max), dims = (2, 3))) - is_auto ? cor[lx:end] : cor -end - -r_autocov(x::AbstractVector) = r_crosscov(x, x, true) - -function r_crosscov(x::AbstractVector, y::AbstractVector, is_auto=false) - cv(x,y) = _ccf(is_auto, x, y) - cc_re = cv(real.(x), real.(y)) + cv(imag.(x), imag.(y)) - cc_im = cv(real.(x), imag.(y)) - cv(imag.(x), real.(y)) - cc_re + im*cc_im -end - -r_autocor(x::AbstractVector) = r_crosscor(x, x, true) - -function r_crosscor(x::AbstractVector, y::AbstractVector, is_auto=false) - cc = r_crosscov(x, y, is_auto) - sc(z) = _r_ccf(real.(z), real.(z), 0) + _r_ccf(imag.(z), imag.(z), 0) - cc/sqrt(sc(x)*sc(y)) -end diff --git a/test/signalcorr.jl b/test/signalcorr.jl index 566a8280b..22da15f3e 100644 --- a/test/signalcorr.jl +++ b/test/signalcorr.jl @@ -29,6 +29,8 @@ cmplx_x2 = convert(AbstractVector{Complex}, x2) @test autocov([1:5;]) ≈ [2.0, 0.8, -0.2, -0.8, -0.8] @test autocor([1, 2, 3, 4, 5]) ≈ [1.0, 0.4, -0.1, -0.4, -0.4] +@test_throws MethodError autocov([1 missing 2 3 4 5]) +@test_throws MethodError autocor([1 missing 2 3 4 5]) acovx1 = [0.755284179631112 + 0.0im, -0.005333112170365584 - 0.18633291805458002im, @@ -144,7 +146,7 @@ function yulewalker_qr(v::AbstractVector) x = -A\b end function toeplitz(v::AbstractVector{T}) where T - N=length(v) + N = length(v) A = zeros(T, N - 1, N - 1) for n in 1:N-1 A[n, n+1:end] = conj(v[2:N-n])