diff --git a/docs/src/scalarstats.md b/docs/src/scalarstats.md
index 6348b75f..1c50795b 100644
--- a/docs/src/scalarstats.md
+++ b/docs/src/scalarstats.md
@@ -7,7 +7,6 @@ The package implements functions for computing various statistics over an array
 ```@docs
 mean
 mean!
-middle
 geomean
 harmmean
 genmean
@@ -58,6 +57,7 @@ quantile
 quantile!
 median
 median!
+middle
 ```
 
 ## Mode and Modes
diff --git a/src/Statistics.jl b/src/Statistics.jl
index 28394ed0..fcabdae2 100644
--- a/src/Statistics.jl
+++ b/src/Statistics.jl
@@ -39,6 +39,7 @@ include("moments.jl")
 include("scalarstats.jl")
 include("cov.jl")
 include("partialcor.jl")
+include("toeplitzsolvers.jl")
 include("signalcorr.jl")
 
 ##### mean #####
@@ -727,8 +728,10 @@ function _getnobs(x::AbstractVecOrMat, y::AbstractVecOrMat, vardim::Int)
     return n
 end
 
-_vmean(x::AbstractVector, vardim::Int) = mean(x)
-_vmean(x::AbstractMatrix, vardim::Int) = mean(x, dims=vardim)
+_vmean(x::AbstractVector, vardim::Int, w::Union{AbstractWeights, Nothing}=nothing) =
+    mean(x, weights=w)
+_vmean(x::AbstractMatrix, vardim::Int, w::Union{AbstractWeights, Nothing}=nothing) =
+    mean(x, dims=vardim, weights=w)
 
 # core functions
 
@@ -771,7 +774,7 @@ end
 ## which can't be handled by broadcast
 covm(x::AbstractVector, xmean; corrected::Bool=true) =
     covzm(map(t -> t - xmean, x); corrected=corrected)
-covm(x::AbstractMatrix, xmean, vardim::Int=1; corrected::Bool=true) =
+covm(x::AbstractMatrix, xmean, weights::Nothing=nothing, vardim::Int=1; corrected::Bool=true) =
     covzm(x .- xmean, vardim; corrected=corrected)
 covm(x::AbstractVector, xmean, y::AbstractVector, ymean; corrected::Bool=true) =
     covzm(map(t -> t - xmean, x), map(t -> t - ymean, y); corrected=corrected)
@@ -788,14 +791,24 @@ is scaled with `n-1`, whereas the sum is scaled with `n` if `corrected` is `fals
 cov(x::AbstractVector; corrected::Bool=true) = covm(x, mean(x); corrected=corrected)
 
 """
-    cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true)
+    cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true[, weights::AbstractWeights])
 
 Compute the covariance matrix of the matrix `X` along the dimension `dims`. If `corrected`
 is `true` (the default) then the sum is scaled with `n-1`, whereas the sum is scaled with `n`
 if `corrected` is `false` where `n = size(X, dims)`.
+
+If `weights` is provided, the biased covariance matrix (`corrected=false`)
+is computed by multiplying `scattermat(X, w)` by
+``\\frac{1}{\\sum{w}}`` to normalize. However, the unbiased covariance matrix
+(`corrected=true`) is dependent on the type of weights used:
+* `AnalyticWeights`: ``\\frac{1}{\\sum w - \\sum {w^2} / \\sum w}``
+* `FrequencyWeights`: ``\\frac{1}{\\sum{w} - 1}``
+* `ProbabilityWeights`: ``\\frac{n}{(n - 1) \\sum w}`` where ``n`` equals `count(!iszero, w)`
+* `Weights`: `ArgumentError` (bias correction not supported)
 """
-cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true) =
-    covm(X, _vmean(X, dims), dims; corrected=corrected)
+cov(X::AbstractMatrix; dims::Int=1, corrected::Bool=true,
+    weights::Union{AbstractWeights, Nothing}=nothing) =
+    covm(X, _vmean(X, dims, weights), weights, dims; corrected=corrected)
 
 """
     cov(x::AbstractVector, y::AbstractVector; corrected::Bool=true)
@@ -899,7 +912,8 @@ corzm(x::AbstractMatrix, y::AbstractMatrix, vardim::Int=1) =
 
 corm(x::AbstractVector{T}, xmean) where {T} =
     T === Missing ? missing : one(float(nonmissingtype(T)))
-corm(x::AbstractMatrix, xmean, vardim::Int=1) = corzm(x .- xmean, vardim)
+corm(x::AbstractMatrix, xmean, weights::Nothing=nothing, vardim::Int=1) =
+    corzm(x .- xmean, vardim)
 function corm(x::AbstractVector, mx, y::AbstractVector, my)
     require_one_based_indexing(x, y)
     n = length(x)
@@ -936,11 +950,13 @@ cor(x::AbstractVector{T}) where {T} =
     T === Missing ? missing : one(float(nonmissingtype(T)))
 
 """
-    cor(X::AbstractMatrix; dims::Int=1)
+    cor(X::AbstractMatrix; dims::Int=1[, weights::AbstractWeights])
 
 Compute the Pearson correlation matrix of the matrix `X` along the dimension `dims`.
+The weighted correlation is computed if `weights` is provided.
 """
-cor(X::AbstractMatrix; dims::Int=1) = corm(X, _vmean(X, dims), dims)
+cor(X::AbstractMatrix; dims::Int=1, weights::Union{AbstractWeights, Nothing}=nothing) =
+    corm(X, _vmean(X, dims, weights), weights, dims)
 
 """
     cor(x::AbstractVector, y::AbstractVector)
diff --git a/src/cov.jl b/src/cov.jl
index d1e51318..72333e0a 100644
--- a/src/cov.jl
+++ b/src/cov.jl
@@ -32,11 +32,9 @@ _unscaled_covzm(x::DenseMatrix, wv::AbstractWeights, dims::Integer) =
     _symmetrize!(unscaled_covzm(x, _scalevars(x, values(wv), dims), dims))
 
 """
-    scattermat(X, [wv::AbstractWeights]; mean=nothing, dims=1)
+    scattermat(X; mean=nothing, dims=1[, weights::AbstractWeights])
 
 Compute the scatter matrix, which is an unnormalized covariance matrix.
-A weighting vector `wv` can be specified to weight
-the estimate.
 
 # Arguments
 * `mean=nothing`: a known mean value. `nothing` indicates that the mean is
@@ -45,62 +43,33 @@ the estimate.
 * `dims=1`: the dimension along which the variables are organized.
   When `dims = 1`, the variables are considered columns with observations in rows;
   when `dims = 2`, variables are in rows with observations in columns.
-"""
-function scattermat end
-
-
-"""
-    cov(X, w::AbstractWeights, vardim=1; mean=nothing,  corrected=false)
-
-Compute the weighted covariance matrix. Similar to `var` and `std` the biased covariance
-matrix (`corrected=false`) is computed by multiplying `scattermat(X, w)` by
-``\\frac{1}{\\sum{w}}`` to normalize. However, the unbiased covariance matrix
-(`corrected=true`) is dependent on the type of weights used:
-* `AnalyticWeights`: ``\\frac{1}{\\sum w - \\sum {w^2} / \\sum w}``
-* `FrequencyWeights`: ``\\frac{1}{\\sum{w} - 1}``
-* `ProbabilityWeights`: ``\\frac{n}{(n - 1) \\sum w}`` where ``n`` equals `count(!iszero, w)`
-* `Weights`: `ArgumentError` (bias correction not supported)
-"""
-cov
-
-scattermat(x::DenseMatrix; mean=nothing, dims::Int=1) =
-    _scattermatm(x, mean, dims)
-_scattermatm(x::DenseMatrix, ::Nothing, dims::Int) =
-    _unscaled_covzm(x .- mean(x, dims=dims), dims)
-_scattermatm(x::DenseMatrix, mean, dims::Int=1) =
+* `weights`: optional weights for observations.
+"""
+scattermat(x::DenseMatrix; mean=nothing, dims::Int=1,
+           weights::Union{AbstractWeights, Nothing}=nothing) =
+    _scattermatm(x, weights, mean, dims)
+_scattermatm(x::DenseMatrix, weights::Nothing, mean::Nothing, dims::Int) =
+    _unscaled_covzm(x .- Statistics.mean(x, dims=dims), dims)
+_scattermatm(x::DenseMatrix, weights::Nothing, mean, dims::Int=1) =
     _unscaled_covzm(x .- mean, dims)
 
-scattermat(x::DenseMatrix, wv::AbstractWeights; mean=nothing, dims::Int=1) =
-    _scattermatm(x, wv, mean, dims)
-_scattermatm(x::DenseMatrix, wv::AbstractWeights, ::Nothing, dims::Int) =
-    _unscaled_covzm(x .- mean(x, wv, dims=dims), wv, dims)
-_scattermatm(x::DenseMatrix, wv::AbstractWeights, mean, dims::Int) =
-    _unscaled_covzm(x .- mean, wv, dims)
+_scattermatm(x::DenseMatrix, weights::AbstractWeights, mean::Nothing, dims::Int) =
+    _unscaled_covzm(x .- Statistics.mean(x, weights=weights, dims=dims), weights, dims)
+_scattermatm(x::DenseMatrix, weights::AbstractWeights, mean, dims::Int) =
+    _unscaled_covzm(x .- mean, weights, dims)
 
 ## weighted cov
-covm(x::DenseMatrix, mean, w::AbstractWeights, dims::Int=1;
+covm(x::DenseMatrix, mean, weights::AbstractWeights, dims::Int=1;
      corrected::Bool=true) =
-    rmul!(scattermat(x, w, mean=mean, dims=dims), varcorrection(w, depcheck(:covm, corrected)))
-
+    rmul!(scattermat(x, weights=weights, mean=mean, dims=dims),
+          varcorrection(weights, corrected))
 
-cov(x::DenseMatrix, w::AbstractWeights, dims::Int=1; corrected::Bool=true) =
-    covm(x, mean(x, w, dims=dims), w, dims; corrected=depcheck(:cov, corrected))
-
-function corm(x::DenseMatrix, mean, w::AbstractWeights, vardim::Int=1)
-    c = covm(x, mean, w, vardim; corrected=false)
-    s = stdm(x, w, mean, vardim; corrected=false)
+function corm(x::DenseMatrix, mean, weights::AbstractWeights, vardim::Int=1)
+    c = covm(x, mean, weights, vardim; corrected=false)
+    s = std(x, mean=mean, weights=weights, dims=vardim, corrected=false)
     cov2cor!(c, s)
 end
 
-"""
-    cor(X, w::AbstractWeights, dims=1)
-
-Compute the Pearson correlation matrix of `X` along the dimension
-`dims` with a weighting `w` .
-"""
-cor(x::DenseMatrix, w::AbstractWeights, dims::Int=1) =
-    corm(x, mean(x, w, dims=dims), w, dims)
-
 """
     cov2cor(C, s)
 
@@ -156,7 +125,8 @@ cov(ce::CovarianceEstimator, x::AbstractVector, y::AbstractVector) =
     error("cov is not defined for $(typeof(ce)), $(typeof(x)) and $(typeof(y))")
 
 """
-    cov(ce::CovarianceEstimator, X::AbstractMatrix, [w::AbstractWeights]; mean=nothing, dims::Int=1)
+    cov(ce::CovarianceEstimator, X::AbstractMatrix; mean=nothing, dims::Int=1,
+        [weights::AbstractWeights])
 
 Compute the covariance matrix of the matrix `X` along dimension `dims`
 using estimator `ce`. A weighting vector `w` can be specified.
@@ -170,18 +140,16 @@ The keyword argument `mean` can be:
   * when `dims=2`, an `AbstractVector` of length `N` or an `AbstractMatrix`
     of size `(N,1)`.
 """
-cov(ce::CovarianceEstimator, X::AbstractMatrix; mean=nothing, dims::Int=1) =
+cov(ce::CovarianceEstimator, X::AbstractMatrix; mean=nothing, dims::Int=1,
+    weights::Union{AbstractWeights, Nothing}=nothing) =
     error("cov is not defined for $(typeof(ce)) and $(typeof(X))")
 
-cov(ce::CovarianceEstimator, X::AbstractMatrix, w::AbstractWeights; mean=nothing, dims::Int=1) =
-    error("cov is not defined for $(typeof(ce)), $(typeof(X)) and $(typeof(w))")
-
 """
     SimpleCovariance(;corrected::Bool=false)
 
 Simple covariance estimator. Estimation calls `cov(x; corrected=corrected)`,
-`cov(x, y; corrected=corrected)` or `cov(X, w, dims; corrected=corrected)`
-where `x`, `y` are vectors, `X` is a matrix and `w` is a weighting vector.
+`cov(x, y; corrected=corrected)` or `cov(X, dims=dims, weights=weights, corrected=corrected)`
+where `x`, `y` are vectors, `X` is a matrix and `weights` is a weighting vector.
 """
 struct SimpleCovariance <: CovarianceEstimator
     corrected::Bool
@@ -194,20 +162,13 @@ cov(sc::SimpleCovariance, x::AbstractVector) =
 cov(sc::SimpleCovariance, x::AbstractVector, y::AbstractVector) =
     cov(x, y; corrected=sc.corrected)
 
-function cov(sc::SimpleCovariance, X::AbstractMatrix; dims::Int=1, mean=nothing)
-    dims ∈ (1, 2) || throw(ArgumentError("Argument dims can only be 1 or 2 (given: $dims)"))
-    if mean === nothing
-        return cov(X; dims=dims, corrected=sc.corrected)
-    else
-        return covm(X, mean, dims, corrected=sc.corrected)
-    end
-end
-
-function cov(sc::SimpleCovariance, X::AbstractMatrix, w::AbstractWeights; dims::Int=1, mean=nothing)
+function cov(sc::SimpleCovariance, X::AbstractMatrix;
+             dims::Int=1,
+             weights::Union{AbstractWeights, Nothing}=nothing,
+             mean=nothing)
     dims ∈ (1, 2) || throw(ArgumentError("Argument dims can only be 1 or 2 (given: $dims)"))
     if mean === nothing
-        return cov(X, w, dims, corrected=sc.corrected)
-    else
-        return covm(X, mean, w, dims, corrected=sc.corrected)
+        mean = Statistics.mean(X, dims=dims, weights=weights)
     end
+    return covm(X, mean, weights, dims, corrected=sc.corrected)
 end
diff --git a/test/cov.jl b/test/cov.jl
index ab310276..b41fe5ce 100644
--- a/test/cov.jl
+++ b/test/cov.jl
@@ -1,9 +1,9 @@
-using StatsBase
+using Statistics
 using LinearAlgebra, Random, Test
 
 struct EmptyCovarianceEstimator <: CovarianceEstimator end
 
-@testset "StatsBase.Covariance" begin
+@testset "Covariance" begin
 weight_funcs = (weights, aweights, fweights, pweights)
 
 @testset "$f" for f in weight_funcs
@@ -24,8 +24,8 @@ weight_funcs = (weights, aweights, fweights, pweights)
     wv1 = f(w1)
     wv2 = f(w2)
 
-    Z1w = X .- mean(X, wv1, dims=1)
-    Z2w = X .- mean(X, wv2, dims=2)
+    Z1w = X .- mean(X, weights=wv1, dims=1)
+    Z2w = X .- mean(X, weights=wv2, dims=2)
 
     ## reference results
 
@@ -45,79 +45,44 @@ weight_funcs = (weights, aweights, fweights, pweights)
         @test scattermat(X)         ≈ S1
         @test scattermat(X, dims=2) ≈ S2
 
-        @test StatsBase.scattermat(X, mean=0)         ≈ Sz1
-        @test StatsBase.scattermat(X, mean=0, dims=2) ≈ Sz2
+        @test scattermat(X, mean=0)         ≈ Sz1
+        @test scattermat(X, mean=0, dims=2) ≈ Sz2
 
-        @test StatsBase.scattermat(X, mean=mean(X, dims=1))         ≈ S1
-        @test StatsBase.scattermat(X, mean=mean(X, dims=2), dims=2) ≈ S2
+        @test scattermat(X, mean=mean(X, dims=1))         ≈ S1
+        @test scattermat(X, mean=mean(X, dims=2), dims=2) ≈ S2
 
-        @test StatsBase.scattermat(X, mean=zeros(1,8))       ≈ Sz1
-        @test StatsBase.scattermat(X, mean=zeros(3), dims=2) ≈ Sz2
+        @test scattermat(X, mean=zeros(1,8))       ≈ Sz1
+        @test scattermat(X, mean=zeros(3), dims=2) ≈ Sz2
 
         @testset "Weighted" begin
-            @test scattermat(X, wv1)         ≈ S1w
-            @test scattermat(X, wv2, dims=2) ≈ S2w
+            @test scattermat(X, weights=wv1)         ≈ S1w
+            @test scattermat(X, weights=wv2, dims=2) ≈ S2w
 
-            @test StatsBase.scattermat(X, wv1, mean=0)         ≈ Sz1w
-            @test StatsBase.scattermat(X, wv2, mean=0, dims=2) ≈ Sz2w
+            @test scattermat(X, weights=wv1, mean=0)         ≈ Sz1w
+            @test scattermat(X, weights=wv2, mean=0, dims=2) ≈ Sz2w
 
-            @test StatsBase.scattermat(X, wv1, mean=mean(X, wv1, dims=1))         ≈ S1w
-            @test StatsBase.scattermat(X, wv2, mean=mean(X, wv2, dims=2), dims=2) ≈ S2w
+            @test scattermat(X, weights=wv1, mean=mean(X, weights=wv1, dims=1))         ≈ S1w
+            @test scattermat(X, weights=wv2, mean=mean(X, weights=wv2, dims=2), dims=2) ≈ S2w
 
-            @test StatsBase.scattermat(X, wv1, mean=zeros(1,8))       ≈ Sz1w
-            @test StatsBase.scattermat(X, wv2, mean=zeros(3), dims=2) ≈ Sz2w
+            @test scattermat(X, weights=wv1, mean=zeros(1,8))       ≈ Sz1w
+            @test scattermat(X, weights=wv2, mean=zeros(3), dims=2) ≈ Sz2w
         end
     end
 
     @testset "Uncorrected" begin
         @testset "Weighted Covariance" begin
-            @test cov(X, wv1; corrected=false)    ≈ S1w ./ sum(wv1)
-            @test cov(X, wv2, 2; corrected=false) ≈ S2w ./ sum(wv2)
-
-            @test StatsBase.covm(X, 0, wv1, 1; corrected=false) ≈ Sz1w ./ sum(wv1)
-            @test StatsBase.covm(X, 0, wv2, 2; corrected=false) ≈ Sz2w ./ sum(wv2)
-
-            @test StatsBase.covm(X, mean(X, wv1, dims=1), wv1, 1; corrected=false) ≈ S1w ./ sum(wv1)
-            @test StatsBase.covm(X, mean(X, wv2, dims=2), wv2, 2; corrected=false) ≈ S2w ./ sum(wv2)
-
-            @test StatsBase.covm(X, zeros(1,8), wv1, 1; corrected=false) ≈ Sz1w ./ sum(wv1)
-            @test StatsBase.covm(X, zeros(3), wv2, 2; corrected=false)   ≈ Sz2w ./ sum(wv2)
-        end
-
-        @testset "Mean and covariance" begin
-            (m, C) = mean_and_cov(X; corrected=false)
-            @test m == mean(X, dims=1)
-            @test C == cov(X, dims=1, corrected=false)
-
-            (m, C) = mean_and_cov(X, 1; corrected=false)
-            @test m == mean(X, dims=1)
-            @test C == cov(X, dims=1, corrected = false)
-
-            (m, C) = mean_and_cov(X, 2; corrected=false)
-            @test m == mean(X, dims=2)
-            @test C == cov(X, dims=2, corrected = false)
-
-            (m, C) = mean_and_cov(X, wv1; corrected=false)
-            @test m == mean(X, wv1, dims=1)
-            @test C == cov(X, wv1, 1, corrected=false)
-
-            (m, C) = mean_and_cov(X, wv1, 1; corrected=false)
-            @test m == mean(X, wv1, dims=1)
-            @test C == cov(X, wv1, 1, corrected=false)
-
-            (m, C) = mean_and_cov(X, wv2, 2; corrected=false)
-            @test m == mean(X, wv2, dims=2)
-            @test C == cov(X, wv2, 2, corrected=false)
+            @test cov(X, weights=wv1; corrected=false)    ≈ S1w ./ sum(wv1)
+            @test cov(X, weights=wv2, dims=2; corrected=false) ≈ S2w ./ sum(wv2)
         end
         @testset "Conversions" begin
-            std1 = std(X, wv1, 1; corrected=false)
-            std2 = std(X, wv2, 2; corrected=false)
+            std1 = std(X, weights=wv1, dims=1; corrected=false)
+            std2 = std(X, weights=wv2, dims=2; corrected=false)
 
-            cov1 = cov(X, wv1, 1; corrected=false)
-            cov2 = cov(X, wv2, 2; corrected=false)
+            cov1 = cov(X, weights=wv1, dims=1; corrected=false)
+            cov2 = cov(X, weights=wv2, dims=2; corrected=false)
 
-            cor1 = cor(X, wv1, 1)
-            cor2 = cor(X, wv2, 2)
+            cor1 = cor(X, weights=wv1, dims=1)
+            cor2 = cor(X, weights=wv2, dims=2)
 
             @testset "cov2cor" begin
                 @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
@@ -137,63 +102,25 @@ weight_funcs = (weights, aweights, fweights, pweights)
     @testset "Corrected" begin
         @testset "Weighted Covariance" begin
             if isa(wv1, Weights)
-                @test_throws ArgumentError cov(X, wv1; corrected=true)
-            else
-                var_corr1 = StatsBase.varcorrection(wv1, true)
-                var_corr2 = StatsBase.varcorrection(wv2, true)
-
-                @test cov(X, wv1; corrected=true)    ≈ S1w .* var_corr1
-                @test cov(X, wv2, 2; corrected=true) ≈ S2w .* var_corr2
-
-                @test StatsBase.covm(X, 0, wv1, 1; corrected=true) ≈ Sz1w .* var_corr1
-                @test StatsBase.covm(X, 0, wv2, 2; corrected=true) ≈ Sz2w .* var_corr2
-
-                @test StatsBase.covm(X, mean(X, wv1, dims=1), wv1, 1; corrected=true) ≈ S1w .* var_corr1
-                @test StatsBase.covm(X, mean(X, wv2, dims=2), wv2, 2; corrected=true) ≈ S2w .* var_corr2
-
-                @test StatsBase.covm(X, zeros(1,8), wv1, 1; corrected=true) ≈ Sz1w .* var_corr1
-                @test StatsBase.covm(X, zeros(3), wv2, 2; corrected=true)   ≈ Sz2w .* var_corr2
-            end
-        end
-        @testset "Mean and covariance" begin
-            (m, C) = mean_and_cov(X; corrected=true)
-            @test m == mean(X, dims=1)
-            @test C == cov(X, dims=1, corrected = true)
-
-            (m, C) = mean_and_cov(X, 1; corrected=true)
-            @test m == mean(X, dims=1)
-            @test C == cov(X, dims=1, corrected = true)
-
-            (m, C) = mean_and_cov(X, 2; corrected=true)
-            @test m == mean(X, dims=2)
-            @test C == cov(X, dims=2, corrected = true)
-
-            if isa(wv1, Weights)
-                @test_throws ArgumentError mean_and_cov(X, wv1; corrected=true)
+                @test_throws ArgumentError cov(X, weights=wv1, corrected=true)
             else
-                (m, C) = mean_and_cov(X, wv1; corrected=true)
-                @test m == mean(X, wv1, dims=1)
-                @test C == cov(X, wv1, 1; corrected=true)
-
-                (m, C) = mean_and_cov(X, wv1, 1; corrected=true)
-                @test m == mean(X, wv1, dims=1)
-                @test C == cov(X, wv1, 1; corrected=true)
+                var_corr1 = Statistics.varcorrection(wv1, true)
+                var_corr2 = Statistics.varcorrection(wv2, true)
 
-                (m, C) = mean_and_cov(X, wv2, 2; corrected=true)
-                @test m == mean(X, wv2, dims=2)
-                @test C == cov(X, wv2, 2; corrected=true)
+                @test cov(X, weights=wv1, corrected=true)    ≈ S1w .* var_corr1
+                @test cov(X, weights=wv2, dims=2, corrected=true) ≈ S2w .* var_corr2
             end
         end
         @testset "Conversions" begin
             if !isa(wv1, Weights)
-                std1 = std(X, wv1, 1; corrected=true)
-                std2 = std(X, wv2, 2; corrected=true)
+                std1 = std(X, weights=wv1, dims=1; corrected=true)
+                std2 = std(X, weights=wv2, dims=2; corrected=true)
 
-                cov1 = cov(X, wv1, 1; corrected=true)
-                cov2 = cov(X, wv2, 2; corrected=true)
+                cov1 = cov(X, weights=wv1, dims=1; corrected=true)
+                cov2 = cov(X, weights=wv2, dims=2; corrected=true)
 
-                cor1 = cor(X, wv1, 1)
-                cor2 = cor(X, wv2, 2)
+                cor1 = cor(X, weights=wv1, dims=1)
+                cor2 = cor(X, weights=wv2, dims=2)
 
                 @testset "cov2cor" begin
                     @test cov2cor(cov(X, dims = 1), std(X, dims = 1)) ≈ cor(X, dims = 1)
@@ -205,12 +132,12 @@ weight_funcs = (weights, aweights, fweights, pweights)
                 @testset "cov2cor!" begin
                     tmp_cov1 = copy(cov1)
                     @test !(tmp_cov1 ≈ cor1)
-                    StatsBase.cov2cor!(tmp_cov1, std1)
+                    Statistics.cov2cor!(tmp_cov1, std1)
                     @test tmp_cov1 ≈ cor1
 
                     tmp_cov2 = copy(cov2)
                     @test !(tmp_cov2 ≈ cor2)
-                    StatsBase.cov2cor!(tmp_cov2, std2)
+                    Statistics.cov2cor!(tmp_cov2, std2)
                     @test tmp_cov2 ≈ cor2
                 end
 
@@ -224,12 +151,12 @@ weight_funcs = (weights, aweights, fweights, pweights)
                 @testset "cor2cov!" begin
                     tmp_cor1 = copy(cor1)
                     @test !(tmp_cor1 ≈ cov1)
-                    StatsBase.cor2cov!(tmp_cor1, std1)
+                    Statistics.cor2cov!(tmp_cor1, std1)
                     @test tmp_cor1 ≈ cov1
 
                     tmp_cor2 = copy(cor2)
                     @test !(tmp_cor2 ≈ cov2)
-                    StatsBase.cor2cov!(tmp_cor2, std2)
+                    Statistics.cor2cov!(tmp_cor2, std2)
                     @test tmp_cor2 ≈ cov2
                 end
             end
@@ -237,18 +164,18 @@ weight_funcs = (weights, aweights, fweights, pweights)
     end
 
     @testset "Correlation" begin
-        @test cor(X, f(ones(3)), 1) ≈ cor(X, dims = 1)
-        @test cor(X, f(ones(8)), 2) ≈ cor(X, dims = 2)
-
-        cov1 = cov(X, wv1, 1; corrected=false)
-        std1 = std(X, wv1, 1; corrected=false)
-        cov2 = cov(X, wv2, 2; corrected=false)
-        std2 = std(X, wv2, 2; corrected=false)
-        expected_cor1 = StatsBase.cov2cor!(cov1, std1)
-        expected_cor2 = StatsBase.cov2cor!(cov2, std2)
-
-        @test cor(X, wv1, 1) ≈ expected_cor1
-        @test cor(X, wv2, 2) ≈ expected_cor2
+        @test cor(X, weights=f(ones(3)), dims=1) ≈ cor(X, dims = 1)
+        @test cor(X, weights=f(ones(8)), dims=2) ≈ cor(X, dims = 2)
+
+        cov1 = cov(X, weights=wv1, dims=1, corrected=false)
+        std1 = std(X, weights=wv1, dims=1, corrected=false)
+        cov2 = cov(X, weights=wv2, dims=2, corrected=false)
+        std2 = std(X, weights=wv2, dims=2, corrected=false)
+        expected_cor1 = Statistics.cov2cor!(cov1, std1)
+        expected_cor2 = Statistics.cov2cor!(cov2, std2)
+
+        @test cor(X, weights=wv1, dims=1) ≈ expected_cor1
+        @test cor(X, weights=wv2, dims=2) ≈ expected_cor2
     end
 
     @testset "Abstract covariance estimation" begin
@@ -258,15 +185,19 @@ weight_funcs = (weights, aweights, fweights, pweights)
         for corrected ∈ (false, true)
             scc = SimpleCovariance(corrected=corrected)
             @test_throws ArgumentError cov(scc, X, dims=0)
-            @test_throws ArgumentError cov(scc, X, wv1, dims=0)
+            @test_throws ArgumentError cov(scc, X, weights=wv1, dims=0)
             @test cov(scc, X) ≈ cov(X, corrected=corrected)
-            @test cov(scc, X, mean=Xm1) ≈ StatsBase.covm(X, Xm1, corrected=corrected)
-            @test cov(scc, X, mean=Xm2, dims=2) ≈ StatsBase.covm(X, Xm2, 2, corrected=corrected)
+            @test cov(scc, X, mean=Xm1) ≈ Statistics.covm(X, Xm1, nothing, corrected=corrected)
+            @test cov(scc, X, mean=Xm2, dims=2) ≈ Statistics.covm(X, Xm2, nothing, 2, corrected=corrected)
             if f !== weights || corrected === false
-                @test cov(scc, X, wv1, dims=1) ≈ cov(X, wv1, 1, corrected=corrected)
-                @test cov(scc, X, wv2, dims=2) ≈ cov(X, wv2, 2, corrected=corrected)
-                @test cov(scc, X, wv1, mean=Xm1) ≈ StatsBase.covm(X, Xm1, wv1, corrected=corrected)
-                @test cov(scc, X, wv2, mean=Xm2, dims=2) ≈ StatsBase.covm(X, Xm2, wv2, 2, corrected=corrected)
+                @test cov(scc, X, weights=wv1, dims=1) ≈
+                    cov(X, weights=wv1, dims=1, corrected=corrected)
+                @test cov(scc, X, weights=wv2, dims=2) ≈
+                    cov(X, weights=wv2, dims=2, corrected=corrected)
+                @test cov(scc, X, weights=wv1, mean=Xm1) ≈
+                    Statistics.covm(X, Xm1, wv1, corrected=corrected)
+                @test cov(scc, X, weights=wv2, mean=Xm2, dims=2) ≈
+                    Statistics.covm(X, Xm2, wv2, 2, corrected=corrected)
             end
         end
     end
@@ -276,13 +207,13 @@ end
     est = EmptyCovarianceEstimator()
     wv = fweights(rand(2))
     @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0])
-    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], wv)
+    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], weights=wv)
     @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], dims = 2)
-    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], wv, dims = 2)
+    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], weights=wv, dims = 2)
     @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], mean = nothing)
-    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], wv, mean = nothing)
+    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], weights=wv, mean = nothing)
     @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], dims = 2, mean = nothing)
-    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], wv, dims = 2, mean = nothing)
+    @test_throws ErrorException cov(est, [1.0 2.0; 3.0 4.0], weights=wv, dims = 2, mean = nothing)
     @test_throws ErrorException cov(est, [1.0, 2.0], [3.0, 4.0])
     @test_throws ErrorException cov(est, [1.0, 2.0])
 
@@ -296,4 +227,4 @@ end
         @test cov(scc, x, y) ≈ cov(x, y; corrected=corrected)
     end
 end
-end # @testset "StatsBase.Covariance"
+end # @testset "Covariance"
diff --git a/test/partialcor.jl b/test/partialcor.jl
index 77ae3cba..b23458b9 100644
--- a/test/partialcor.jl
+++ b/test/partialcor.jl
@@ -1,4 +1,4 @@
-using StatsBase
+using Statistics
 using Test
 
 wechsler = Float32[
diff --git a/test/runtests.jl b/test/runtests.jl
index 62a12603..e40704c6 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -385,7 +385,7 @@ Y = [6.0  2.0;
         @inferred cov(x1, corrected=cr)
 
         @test cov(X) == Statistics.covm(X, mean(X, dims=1))
-        C = zm ? Statistics.covm(X, 0, vd, corrected=cr) :
+        C = zm ? Statistics.covm(X, 0, nothing, vd, corrected=cr) :
                  cov(X, dims=vd, corrected=cr)
         @test size(C) == (k, k)
         @test C ≈ Cxx
@@ -471,7 +471,7 @@ end
         @inferred cor(x1)
 
         @test cor(X) == Statistics.corm(X, mean(X, dims=1))
-        C = zm ? Statistics.corm(X, 0, vd) : cor(X, dims=vd)
+        C = zm ? Statistics.corm(X, 0, nothing, vd) : cor(X, dims=vd)
         @test size(C) == (k, k)
         @test C ≈ Cxx
         @inferred cor(X, dims=vd)
diff --git a/test/signalcorr.jl b/test/signalcorr.jl
index c8ee35dd..78a14f4f 100644
--- a/test/signalcorr.jl
+++ b/test/signalcorr.jl
@@ -4,7 +4,7 @@
 # The reference results are generated from R.
 #
 
-using StatsBase
+using Statistics
 using Test
 
 # random data for testing