From f1d6e89b4ef6f24f5cb710c4f71d35afabf76af0 Mon Sep 17 00:00:00 2001 From: Lukas Date: Thu, 6 Apr 2023 00:10:02 +0200 Subject: [PATCH] Add support for creating correlated variables from a covariance or correlation matrix (#137) * Add functions to calculate the correlation and covariance matrix * Add functions to create correlated measurements from covariance matrix * Remove TODO entry --- Project.toml | 1 + benchmark/Manifest.toml | 2 +- docs/Project.toml | 2 + docs/src/examples.md | 40 ++++++++++++++++ docs/src/index.md | 2 +- docs/src/todo.md | 6 --- docs/src/usage.md | 26 +++++++++++ src/Measurements.jl | 5 ++ src/utils.jl | 89 +++++++++++++++++++++++++++++++++++ test/runtests.jl | 100 ++++++++++++++++++++++++++++++++++++++++ 10 files changed, 265 insertions(+), 8 deletions(-) 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]