diff --git a/Project.toml b/Project.toml
index 1d165270..6f491d55 100644
--- a/Project.toml
+++ b/Project.toml
@@ -9,6 +9,7 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 RecipesBase = "3cdcf5f2-1ef4-517c-9805-6587b60abb01"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
 Calculus = "0.4.1, 0.5"
diff --git a/benchmark/Manifest.toml b/benchmark/Manifest.toml
index 34530c07..067de4d1 100644
--- a/benchmark/Manifest.toml
+++ b/benchmark/Manifest.toml
@@ -224,7 +224,7 @@ uuid = "c8ffd9c3-330d-5841-b78e-0817d7145fa1"
 version = "2.28.0+0"
 
 [[deps.Measurements]]
-deps = ["Calculus", "LinearAlgebra", "Printf", "RecipesBase", "Requires"]
+deps = ["Calculus", "LinearAlgebra", "Printf", "RecipesBase", "Requires", "Statistics"]
 path = ".."
 uuid = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 version = "2.9.0"
diff --git a/docs/Project.toml b/docs/Project.toml
index 770a9a87..78f8e3d8 100644
--- a/docs/Project.toml
+++ b/docs/Project.toml
@@ -3,8 +3,10 @@ AutoGrad = "6710c13c-97f1-543f-91c5-74e8f7d95b35"
 Calculus = "49dc2e85-a5d0-5ad3-a950-438e2897f1b9"
 Cuba = "8a292aeb-7a57-582c-b821-06e4c11590b1"
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+LsqFit = "2fda8390-95c7-5789-9bda-21331edee243"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [compat]
diff --git a/docs/src/examples.md b/docs/src/examples.md
index 911d8d06..8b8af0f5 100644
--- a/docs/src/examples.md
+++ b/docs/src/examples.md
@@ -446,6 +446,46 @@ Measurements.value.([complex(87.3 ± 2.9, 64.3 ± 3.0), complex(55.1 ± 2.8, -19
 Measurements.uncertainty.([complex(87.3 ± 2.9, 64.3 ± 3.0), complex(55.1 ± 2.8, -19.1 ± 4.6)])
 ```
 
+Calculating the Covariance and Correlation Matrices
+---------------------------------------------------
+
+Calculate the covariance and correlation matrices of multiple `Measurement`s
+with the functions [`cov(::AbstractVector{<:Measurement})`](@ref) and
+[`cor(::AbstractVector{<:Measurement})`](@ref):
+
+```@repl
+using Measurements
+
+x = measurement(1.0, 0.1)
+y = -2x + 10
+z = -3x
+cov([x, y, z])
+cor([x, y, z])
+```
+
+Creating Correlated Measurements from their Nominal Values and a Covariance Matrix
+----------------------------------------------------------------------------------
+
+Using [`Measurements.correlated_values`](@ref), you can correlate
+the uncertainty of fit parameters given a covariance matrix from
+the fitted results. An example using
+[LsqFit.jl](https://julianlsolvers.github.io/LsqFit.jl/latest/):
+
+```@repl
+using Measurements, LsqFit
+
+@. model(x, p) = p[1] + x * p[2]
+
+xdata = range(0, stop=10, length=20)
+ydata = model(xdata, [1.0 2.0]) .+ 0.01 .* randn.()
+p0 = [0.5, 0.5]
+
+fit = curve_fit(model, xdata, ydata, p0)
+
+coeficients = Measurements.correlated_values(coef(fit), estimate_covar(fit))
+f(x) = model(x, coeficients)
+```
+
 Interplay with Third-Party Packages
 -----------------------------------
 
diff --git a/docs/src/index.md b/docs/src/index.md
index 3385307f..3feb6974 100644
--- a/docs/src/index.md
+++ b/docs/src/index.md
@@ -84,7 +84,7 @@ convenience, a BibTeX entry is provided in the
 file.
 
 Other features are expected to come in the future, see the [How Can I
-Help?](@ref) section and the [TODO](@ref) list.
+Help?](@ref) section.
 
 The `Measurements.jl` package is licensed under the MIT "Expat" License. The
 original author is Mosè Giordano.
diff --git a/docs/src/todo.md b/docs/src/todo.md
index 9cca922d..39163877 100644
--- a/docs/src/todo.md
+++ b/docs/src/todo.md
@@ -5,12 +5,6 @@ The package is developed at
 https://github.com/JuliaPhysics/Measurements.jl. There you can submit bug
 reports, make suggestions, and propose pull requests.
 
-TODO
-----
-
-* Allow defining quantities related by a correlation matrix and correctly
-  propagate uncertainty in this case
-
 How Can I Help?
 ---------------
 
diff --git a/docs/src/usage.md b/docs/src/usage.md
index a14bd10d..a26f5631 100644
--- a/docs/src/usage.md
+++ b/docs/src/usage.md
@@ -206,6 +206,32 @@ value and the uncertainty of `x`, be it a single measurement or an array of
 measurements. They are particularly useful in the case of complex measurements
 or arrays of measurements.
 
+Calculating the Covariance and Correlation Matrices
+---------------------------------------------------
+
+```@docs
+Measurements.Statistics.cov(::AbstractVector{<:Measurement})
+Measurements.Statistics.cor(::AbstractVector{<:Measurement})
+```
+
+The [covariance matrix](https://en.wikipedia.org/wiki/Covariance_matrix) of
+multiple measurements can be calculated using the `cov` function by supplying
+a vector of `Measurement`s. Likewise, the [correlation
+matrix](https://en.wikipedia.org/wiki/Correlation#Correlation_matrices) can be
+calculated using the `cor` function with the same signature.
+
+Creating Correlated Measurements from their Nominal Values and a Covariance Matrix
+----------------------------------------------------------------------------------
+
+```@docs
+Measurements.correlated_values
+```
+
+Given some nominal values with an associated covariance matrix, you can
+construct measurements with a correlated uncertainty. Providing both an
+`AbstractVector{<:Real}` of nominal values and a covariance matrix of type
+`AbstractMatrix{<:Real}`.
+
 Error Propagation of Numbers with Units
 ---------------------------------------
 
diff --git a/src/Measurements.jl b/src/Measurements.jl
index e330192b..b8be0909 100644
--- a/src/Measurements.jl
+++ b/src/Measurements.jl
@@ -23,11 +23,16 @@ module Measurements
 
 # Calculus is used to calculate numerical derivatives in "@uncertain" macro.
 using Calculus
+using Statistics
+
+import Statistics: cor, cov
 
 using Requires
 
 # Functions provided by this package and exposed to users
 export Measurement, measurement, ±
+# Re-export from Statistics
+export cor, cov
 
 # Define the "Derivatives" type, used inside "Measurement" type.  This should be
 # a lightweight and immutable dictionary.
diff --git a/src/utils.jl b/src/utils.jl
index 1043703f..e823111f 100644
--- a/src/utils.jl
+++ b/src/utils.jl
@@ -105,3 +105,92 @@ function uncertainty_components(x::Measurement{T}) where {T<:AbstractFloat}
     end
     return out
 end
+
+"""
+    cov(x::AbstractVector{<:Measurement})
+
+Returns the covariance matrix of a vector of correlated `Measurement`s.
+"""
+function Statistics.cov(x::AbstractVector{Measurement{T}}) where T
+    S = length(x)
+    covariance_matrix = zeros(T, (S, S))
+
+    for (ii, i) = enumerate(eachindex(x)), (jj, j) = Iterators.take(enumerate(eachindex(x)), ii)
+        overlap = keys(x[i].der) ∩ keys(x[j].der)
+        covariance_matrix[ii, jj] = isempty(overlap) ? 0.0 : sum(overlap) do var
+            x[i].der[var] * x[j].der[var] * var[2]^2
+        end
+    end
+
+    return Symmetric(covariance_matrix, :L)
+end
+
+"""
+    cor(x::AbstractVector{<:Measurement})
+
+Returns the correlation matrix of a vector of correlated `Measurement`s.
+"""
+function Statistics.cor(x::AbstractVector{<:Measurement})
+    covariance_matrix = cov(x)
+    standard_deviations = sqrt.(diag(covariance_matrix))
+    return covariance_matrix ./ standard_deviations ./ standard_deviations'
+end
+
+"""
+    Measurements.correlated_values(
+        nominal_values::AbstractVector{<:Real},
+        covariance_matrix::AbstractMatrix{<:Real}
+    )
+
+Returns correlated `Measurement`s given their nominal values
+and their covariance matrix.
+"""
+function correlated_values(
+    nominal_values::AbstractVector{<:Real},
+    covariance_matrix::AbstractMatrix{<:Real},
+)
+    standard_deviations = sqrt.(diag(covariance_matrix))
+
+    standard_deviations_nonzero = deepcopy(standard_deviations)
+    standard_deviations_nonzero[findall(iszero, standard_deviations_nonzero)] .= 1
+
+    correlation_matrix = covariance_matrix ./ standard_deviations_nonzero ./ standard_deviations_nonzero'
+
+    return correlated_values(nominal_values, standard_deviations, correlation_matrix)
+end
+
+"""
+    Measurements.correlated_values(
+        nominal_values::AbstractVector{<:Real},
+        standard_deviations::AbstractVector{<:Real},
+        correlation_matrix::AbstractMatrix{<:Real}
+    )
+
+Returns correlated `Measurement`s given their nominal values,
+their standard deviations and their correlation matrix.
+"""
+function correlated_values(
+    nominal_values::AbstractVector{<:Real},
+    standard_deviations::AbstractVector{<:Real},
+    correlation_matrix::AbstractMatrix{<:Real},
+)
+    eig = eigen(correlation_matrix)::Eigen{Float64}
+
+    variances = eig.values
+
+    # Inverse is equal to transpose for unitary matrices
+    transform = eig.vectors'
+
+    # Avoid numerical errors
+    variances[findall(x -> x < eps(1.0), variances)] .= 0
+    variables = map(variances) do variance
+        measurement(0, sqrt(variance))
+    end
+
+    transform .*= standard_deviations'
+    transform_columns = [view(transform,:,i) for i in axes(transform, 2)]
+
+    values_funcs = [Measurements.result(n, t, variables) for (n, t) in zip(nominal_values, transform_columns)]
+
+    return values_funcs
+end
diff --git a/test/runtests.jl b/test/runtests.jl
index 5b5126c6..570e5da2 100644
--- a/test/runtests.jl
+++ b/test/runtests.jl
@@ -90,6 +90,106 @@ end
     @test length((x - x + y - y + w - w + c).der) == 0
 end
 
+@testset "Covariance matrix" begin
+    x = measurement(1.0, 0.1)
+    y = -2x + 10
+    z = -3x
+    @test @inferred(cov([x, y, z])) ≈ [0.01 -0.02 -0.03; -0.02 0.04 0.06; -0.03 0.06 0.09]
+
+    u = measurement(1, 0.05)
+    v = measurement(10, 0.1)
+    w = u + 2v
+    @test @inferred(cov([u, v, w])) ≈ [0.0025 0.0 0.0025; 0.0 0.01 0.02; 0.0025 0.02 0.0425]
+end
+
+@testset "Correlation matrix" begin
+    x = measurement(1.0, 0.1)
+    y = -2x + 10
+    z = -3x
+    @test @inferred(cor([x, y, z])) ≈ [1.0 -1.0 -1.0; -1.0 1.0 1.0; -1.0 1.0 1.0]
+
+    u = measurement(1, 0.05)
+    v = measurement(10, 0.1)
+    w = u + 2v
+    @test @inferred(cor([u, v, w])) ≈ [1.0 0.0 0.24253562503633297; 0.0 1.0 0.9701425001453319; 0.24253562503633297 0.9701425001453319 1.0] atol=1e-15
+end
+
+@testset "Correlated values" begin
+    @testset "Simple" begin
+        u_in = measurement(1, 0.1)
+        cov_matrix = @inferred cov([u_in])
+
+        u_out, = @inferred Measurements.correlated_values([1], cov_matrix)
+        @test 2u_in ≈ 2u_out
+    end
+
+    @testset "Covariances" begin
+        x_in = measurement(1, 0.1)
+        y_in = measurement(2, 0.3)
+        z_in = -3x_in + y_in
+
+        xyz_in = [x_in, y_in, z_in]
+
+        cov_matrix = cov([x_in, y_in, z_in])
+
+        @test Measurements.uncertainty.([x_in, y_in, z_in]) .^ 2 ≈ diag(cov_matrix)
+
+        x_out, y_out, z_out = xyz_out = Measurements.correlated_values(Measurements.value.([x_in, y_in, z_in]), cov_matrix)
+
+        @test xyz_out ≈ xyz_in
+        @test z_out ≈ -3x_out + y_out
+        @test 0 ± 0 ≈ -3x_out + y_out - z_out atol=1e-15
+    end
+
+    @testset "Functional relations" begin
+        u_in = measurement(1, 0.05)
+        v_in = measurement(10, 0.1)
+        w_in = u_in + 2v_in
+
+        cov_matrix = cov([u_in, v_in, w_in])
+
+        (u_out, v_out, w_out) = Measurements.correlated_values(Measurements.value.([u_in, v_in, w_in]), cov_matrix)
+
+        @test u_in ≈ u_out
+        @test v_in ≈ v_out
+        @test w_in ≈ w_out
+        @test w_out ≈ u_out + 2v_out
+        @test 0 ± 0 ≈ u_out + 2v_out - w_out atol=2e-8
+
+        corr_matrix = cor([u_out, v_out, w_out])
+        @test corr_matrix[1,1] ≈ 1
+        @test corr_matrix[2,3] ≈ 2Measurements.uncertainty(v_out)/Measurements.uncertainty(w_out)
+    end
+
+    @testset "Numerical robustnes" begin
+        cov_matrix = [1e-70 0.9e-70 -3e-34; 0.9e-70 1e-70 -3e-34; -3e-34 -3e-34 1e10]
+        variables = Measurements.correlated_values([0, 0, 0], cov_matrix)
+
+        @test cov_matrix[1,1] ≈ Measurements.uncertainty(variables[1])^2
+        @test cov_matrix[2,2] ≈ Measurements.uncertainty(variables[2])^2
+        @test cov_matrix[3,3] ≈ Measurements.uncertainty(variables[3])^2
+    end
+
+    @testset "0 variance" begin
+        x_in = measurement(1, 0)
+        y_in = measurement(2, 0)
+        z_in = measurement(3, 5)
+        cov_matrix = cov([x_in, y_in, z_in])
+        nom_values = Measurements.value.([x_in, y_in, z_in])
+        variables = Measurements.correlated_values(nom_values, cov_matrix)
+
+        for (variable, nom_value, variance) in zip(
+            variables, nom_values, diag(cov_matrix))
+
+            @test Measurements.value(variable) ≈ nom_value
+            @test Measurements.uncertainty(variable)^2 ≈ variance
+        end
+        
+        @test cov_matrix ≈ cov(variables)
+        @test [x_in, y_in, z_in] ≈ variables
+    end
+end
+
 @testset "Contributions to uncertainty" begin
     @test sort(collect(values(@inferred(Measurements.uncertainty_components(w * x * y))))) ≈
         [0.2, 0.3, 0.36]