Add functions to calculate the correlation and covariance matrix

docs/src/examples.md

@@ -446,6 +446,21 @@ 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 between a vector of
+measurements using the functions [`Measurements.covariance_matrix`](@ref)
+and [`Measurements.correlation_matrix`](@ref):
+
+```@repl
+x = measurement(1.0, 0.1)
+y = -2x + 10
+z = -3x
+Measurements.covariance_matrix([x, y, z])
+Measurements.correlation_matrix([x, y, z])
+```
+
 Interplay with Third-Party Packages
 -----------------------------------
 

docs/src/usage.md

@@ -206,6 +206,21 @@ 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.correlation_matrix
+Measurements.covariance_matrix
+```
+
+Calculating the [covariance
+matrix](https://en.wikipedia.org/wiki/Covariance_matrix) between any number of
+`Measurement`s can be done by providing a `AbstractVector{<:Measurement}` to
+`Measurements.covariance_matrix`. Equivalently, the [correlation
+matrix](https://en.wikipedia.org/wiki/Correlation#Correlation_matrices) can be
+retrieved using `Measurements.covariance_matrix`.
+
 Error Propagation of Numbers with Units
 ---------------------------------------
 

src/utils.jl

@@ -105,3 +105,34 @@ function uncertainty_components(x::Measurement{T}) where {T<:AbstractFloat}
     end
     return out
 end
+
+"""
+    Measurements.covariance_matrix(n::AbstractVector{<:Measurement})
+
+Returns the covariance matrix of a vector of correlated `Measurement`s.
+"""
+function covariance_matrix(n::AbstractVector{<:Measurement})
+    S = length(n)
+    covariance_matrix = zeros(Float64, (S, S))
+
+    for i = 1:S, j = 1:i
+        overlap = keys(n[i].der) ∩ keys(n[j].der)
+        covariance_matrix[i, j] = isempty(overlap) ? 0.0 : sum(overlap) do var
+            n[i].der[var] * n[j].der[var] * var[2]^2
+        end
+    end
+
+    return Symmetric(covariance_matrix, :L)
+end
+
+"""
+    Measurements.correlation_matrix(n::AbstractVector{<:Measurement})
+
+Returns the correlation matrix of a vector of correlated `Measurement`s.
+"""
+function correlation_matrix(n::AbstractVector{<:Measurement})
+    covariance_matrix = Measurements.covariance_matrix(n)
+    standard_deviations = sqrt.(diag(covariance_matrix))
+    correlation_matrix = covariance_matrix ./ (standard_deviations * standard_deviations')
+    return correlation_matrix
+end

test/runtests.jl

@@ -90,6 +90,30 @@ 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(Measurements.covariance_matrix([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(Measurements.covariance_matrix([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(Measurements.correlation_matrix([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(Measurements.correlation_matrix([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 "Contributions to uncertainty" begin
     @test sort(collect(values(@inferred(Measurements.uncertainty_components(w * x * y))))) ≈
         [0.2, 0.3, 0.36]