Add functions to create correlated measurements from covariance matrix

docs/src/examples.md

@@ -461,6 +461,30 @@ Measurements.covariance_matrix([x, y, z])
 Measurements.correlation_matrix([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/):
+
+```julia
+using LsqFit
+using Measurements
+
+@. 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(length(xdata))
+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
 -----------------------------------
 

docs/src/usage.md

@@ -221,6 +221,18 @@ matrix](https://en.wikipedia.org/wiki/Covariance_matrix) between any number of
 matrix](https://en.wikipedia.org/wiki/Correlation#Correlation_matrices) can be
 retrieved using `Measurements.covariance_matrix`.
 
+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
 ---------------------------------------
 

src/utils.jl

@@ -136,3 +136,62 @@ function correlation_matrix(n::AbstractVector{<:Measurement})
     correlation_matrix = covariance_matrix ./ (standard_deviations * standard_deviations')
     return correlation_matrix
 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},
+)
+    variances = diag(covariance_matrix)
+
+    standard_deviations = sqrt.(variances)
+    standard_deviations[findall(iszero, standard_deviations)] .= 1
+
+    correlation_matrix = covariance_matrix ./ standard_deviations ./ standard_deviations'
+
+    return correlated_values(nominal_values, sqrt.(variances), 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

test/runtests.jl

@@ -114,6 +114,76 @@ end
     @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 "Correlated values" begin
+    @testset "Simple" begin
+        u = measurement(1, 0.1)
+        cov = @inferred Measurements.covariance_matrix([u])
+
+        u2, = @inferred Measurements.correlated_values([1], cov)
+        @test 2u ≈ 2u2
+    end
+
+    @testset "Covariances" begin
+        x = measurement(1, 0.1)
+        y = measurement(2, 0.3)
+        z = -3x + y
+
+        xyz = [x, y, z]
+
+        covs = Measurements.covariance_matrix(xyz)
+
+        @test Measurements.uncertainty.(xyz) .^ 2 ≈ diag(covs)
+
+        x_new, y_new, z_new = xyz_new = Measurements.correlated_values(Measurements.value.(xyz), covs)
+
+        @test xyz_new ≈ xyz
+        @test z_new ≈ -3x_new + y_new
+    end
+
+    @testset "Functional relations" begin
+        u = measurement(1, 0.05)
+        v = measurement(10, 0.1)
+        w = u + 2v
+
+        cov_matrix = Measurements.covariance_matrix([u, v, w])
+
+        (u2, v2, w2) = Measurements.correlated_values(Measurements.value.([u, v, w]), cov_matrix)
+
+        @test u ≈ u2
+        @test v ≈ v2
+        @test w ≈ w2
+        @test 0 ± 0 ≈ w2 - (u2 + 2v2) atol=1e-7
+
+        corr_matrix = Measurements.correlation_matrix([u, v, w])
+        @test corr_matrix[1,1] ≈ 1
+        @test corr_matrix[2,3] ≈ 2Measurements.uncertainty(v)/Measurements.uncertainty(w)
+    end
+
+    @testset "Numerical robustnes" begin
+        cov = [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)
+
+        @test cov[1,1] ≈ Measurements.uncertainty(variables[1])^2
+        @test cov[2,2] ≈ Measurements.uncertainty(variables[2])^2
+        @test cov[3,3] ≈ Measurements.uncertainty(variables[3])^2
+    end
+
+    @testset "0 variance" begin
+        cov = Diagonal([0, 0, 10])
+        nom_values = [1, 2, 3]
+        variables = Measurements.correlated_values(nom_values, cov)
+
+        for (variable, nom_value, variance) in zip(
+            variables, nom_values, diag(cov))
+
+            @test Measurements.value(variable) ≈ nom_value
+            @test Measurements.uncertainty(variable)^2 ≈ variance
+        end
+
+        @test cov ≈ Measurements.covariance_matrix(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]