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