diff --git a/docs/src/examples.md b/docs/src/examples.md index 7af46774..18420286 100644 --- a/docs/src/examples.md +++ b/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 ----------------------------------- diff --git a/docs/src/usage.md b/docs/src/usage.md index 25220608..ce05331c 100644 --- a/docs/src/usage.md +++ b/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 --------------------------------------- diff --git a/src/utils.jl b/src/utils.jl index ad35b14d..2ef1c7f4 100644 --- a/src/utils.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 49b67287..08024056 100644 --- a/test/runtests.jl +++ b/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]