Measurements¶

Correlation Between Variables¶

The fact that two or more measurements are correlated means that there is some sort of relationship beetween them. In the context of measurements and error propagation theory, the term “correlation” is very broad and can indicate different things. Among others, there may be some dependence between uncertainties of different measurements with different values, or a dependence between the values of two measurements while their uncertainties are different.

Here, for correlation we mean the most simple case of functional relationship: if \(x = \bar{x} \pm \sigma_x\) is an independent measurement, a quantity \(y = f(x) = \bar{y} \pm \sigma_y\) that is function of \(x\) is not like an independent measurement but is a quantity that depends on \(x\), so we say that \(y\) is correlated with \(x\). The package Measurements.jl is able to handle this type of correlation when propagating the uncertainty for operations and functions taking two or more arguments. As a result, \(x - x = 0 \pm 0\) and \(x/x = 1 \pm 0\). If this correlation was not accounted for, you would always get non-zero uncertainties even for these operations that have exact results. Two truly different measurements that only by chance share the same nominal value and uncertainty are not treated as correlated.

Propagate Uncertainty for Arbitrary Functions¶

@uncertain f(x, ...)

Existing functions implemented exclusively in Julia that accept AbstractFloat arguments will work out-of-the-box with Measurement objects as long as they internally use functions already supported by this package. However, there are functions that take arguments that are specific subtypes of AbstractFloat, or are implemented in such a way that does not play nicely with Measurement variables.

The package provides the @uncertain macro that overcomes this limitation and further extends the power of Measurements.jl.

This macro allows you to propagate uncertainty in arbitrary functions, including those based on C/Fortran calls, that accept any number of real arguments. The macro exploits derivative and gradient functions from Calculus package in order to perform numerical differentiation.

Derivative and Gradient¶

Measurements.derivative(y::Measurement, x::Measurement)¶

Measurements.gradient(y::Measurement, x::AbstractArray{Measurement})¶

In order to propagate the uncertainties, Measurements.jl keeps track of the partial derivative of an expression with respect to all independent measurements from which the expression comes. For this reason, the package provides two convenient functions, Measurements.derivative and Measurements.gradient, to get the partial derivative and the gradient of an expression with respect to independent measurements.

Standard Score¶

stdscore(measure::Measurement, expected_value::Real) → standard_score¶

The stdscore function is available to calculate the standard score between a measurement and its expected value.

Weighted Average¶

weightedmean(iterable) → weighted_mean¶

weightedmean function gives the weighted mean of a set of measurements using inverses of variances as weights. Use mean for the simple arithmetic mean.

Correlation Between Variables¶

Here you can see examples of how functionally correlated variables are treated within the package:

x = 8.4 ± 0.7
x - x
# => 0.0 ± 0.0
x/x
# => 1.0 ± 0.0
x*x*x - x^3
# => 0.0 ± 0.0
sin(x)/cos(x) - tan(x)
# => -2.220446049250313e-16 ± 0.0
# They are equal within numerical accuracy
y = -5.9 ± 0.2
beta(x, y) - gamma(x)*gamma(y)/gamma(x + y)
# => 0.0 ± 3.979039320256561e-14

You will get similar results for a variable that is a function of an already existing Measurement object:

u = 2x
(x + x) - u
# => 0.0 ± 0.0
u/2x
# => 1.0 ± 0.0
u^3 - 8x^3
# => 0.0 ± 0.0
cos(x)^2 - (1 + cos(u))/2
# => 0.0 ± 0.0

A variable that has the same nominal value and uncertainty as u above but is not functionally correlated with x will give different outcomes:

# Define a new measurement but with same nominal value
# and uncertainty as u, so v is not correlated with x
v = 16.8 ± 1.4
(x + x) - v
# => 0.0 ± 1.979898987322333
v/2x
# => 1.0 ± 0.11785113019775792
v^3 - 8x^3
# => 0.0 ± 1676.4200705455657
cos(x)^2 - (1 + cos(v))/2
# => 0.0 ± 0.8786465354843539

@uncertain Macro¶

Macro @uncertain can be used to propagate uncertainty in arbitrary real or complex functions of real arguments, including functions not natively supported by this package.

@uncertain (x -> complex(zeta(x), exp(eta(x)^2)))(2 ± 0.13)
# => (1.6449340668482273 ± 0.12188127308075564) + (1.9668868646839253 ± 0.042613944993428333)im
@uncertain log(9.4 ± 1.3, 58.8 ± 3.7)
# => 1.8182372640255153 ± 0.11568300475873611
log(9.4 ± 1.3, 58.8 ± 3.7)
# => 1.8182372640255153 ± 0.11568300475593848

You usually do not need to define a wrapping function before using it. In the case where you have to define a function, like in the first line of previous examples, anonymous functions allow you to do it in a very concise way.

The macro works with functions calling C/Fortran functions as well. For example, Cuba.jl package performs numerical integration by wrapping the C Cuba library. You can define a function to numerically compute with Cuba.jl the integral defining the error function and pass it to @uncertain macro. Compare the result with that of the erf function, natively supported in Measurements.jl package

using Cuba
cubaerf(x::Real) =
    2x/sqrt(pi)*Cuhre((t, f) -> f[1] = exp(-abs2(t[1]*x)), 1, 1)[1][1]
@uncertain cubaerf(0.5 ± 0.01)
# => 0.5204998778130466 ± 0.008787825789336267
erf(0.5 ± 0.01)
# => 0.5204998778130465 ± 0.008787825789354449

Also here you can use an anonymous function instead of defining the cubaerf function, do it as an exercise.

atan2(10, 13.5 ± 0.8)
# => 0.6375487981386927 ± 0.028343666961913202
@uncertain atan2(10 ± 0, 13.5 ± 0.8)
# => 0.6375487981386927 ± 0.028343666962347438

@uncertain (x -> airyx(1, x))(4.8 ± 0.2)
# => -0.42300740589773583 ± 0.004083414330362105

Complex Measurements¶

Here are a few examples about uncertainty propagation of complex-valued measurements.

u = complex(32.7 ± 1.1, -3.1 ± 0.2)
v = complex(7.6 ± 0.9, 53.2 ± 3.4)
2u+v
# => (73.0 ± 2.3769728648009427) + (47.0 ± 3.4234485537247377)*im
sqrt(u*v)
# => (33.004702573592 ± 1.0831254428098636) + (25.997507418428984 ± 1.1082833691607152)*im
gamma(u/v)
# => (-0.25050193836584694 ± 0.011473098558745594) + (1.2079738483289788 ± 0.133606565257322)*im

You can also verify the Euler’s formula

cis(u)
# => (6.27781144696534 ± 23.454542573739754) + (21.291738410228678 ± 8.112997844397572)*im
cos(u) + sin(u)*im
# => (6.277811446965339 ± 23.454542573739754) + (21.291738410228678 ± 8.112997844397572)im

Arbitrary Precision Calculations¶

Arbitrary precision calculations involving quantities that are intrinsically imprecise may not be very useful, but Julia natively supports this type of arithmetic and so Measurements.jl does. You only have to create Measurement objects with nominal value and/or uncertainty of type BigFloat (or BigInt as well, actually):

a = BigInt(3) ± 0.01
b = 4 ± 0.03
hypot(a, b)
# => 5.000000000000000000000000000000000000000000000000000000000000000000000000000000 ± 2.473863375370596246756154793364399326509001412701084709723336101627452857843757e-02
log(2a)^b
# => 1.030668097314957384421465902631648727333270687596715387736946157489404400228445e+01 ± 1.959580475953079233643030915452927748488408893913287402297342303952280925878254e-01

Arrays of Measurements¶

You can create arrays of Measurement objects and perform mathematical operations on them in the most natural way possible:

A = [1.03 ± 0.14, 2.88 ± 0.35, 5.46 ± 0.97]
B = [0.92 ± 0.11, 3.14 ± 0.42, 4.67 ± 0.58]
exp(sqrt(B)) - log(A)
# => 3-element Array{Measurements.Measurement{Float64},1}:
#     2.5799612193837493 ± 0.20215123893809778
#     4.824843081566397 ± 0.7076631767039828
#     6.982522998771525 ± 1.178287422979362
cos(A).^2 + sin(A).^2
# 3-element Array{Measurements.Measurement{Float64},1}:
#     1.0 ± 0.0
#     1.0 ± 0.0
#     1.0 ± 0.0
sum(A)
# => 9.370000000000001 ± 1.0406728592598156
mean(B - A)
# => -0.21333333333333326 ± 0.42267665603337445

Derivative and Gradient¶

In order to propagate the uncertainties, Measurements.jl keeps track of the partial derivative of an expression with respect to all independent measurements from which the expression comes. The package provides two convenient functions, Measurements.derivative and Measurements.gradient, that return the partial derivative and the gradient of an expression with respect to independent measurements.

x = 98.1 ± 12.7
y = 105.4 ± 25.6
z = 78.3 ± 14.1
Measurements.derivative(2x - 4y, x)
# => 2.0
Measurements.derivative(2x - 4y, y)
# => -4.0
Measurements.gradient(2x - 4y, [x, y, z])
# => 3-element Array{Float64,1}:
#      2.0
#     -4.0
#      0.0  # The expression does not depend on z

w = y^(3//4)*log(y) + 3x - cos(y/x)
# => 447.0410543780643 ± 52.41813324207829
(Measurements.gradient(w, [x, y]).*[a.err for a in (x,y)]).^2
# => 2-element Array{Any,1}:
#     1442.31
#     1305.36

vecnorm(Measurements.gradient(w, [x, y]).*[a.err for a in (x,y)])
# => 52.41813324207829

stdscore Function¶

You can get the distance in number of standard deviations between a measurement and its expected value (this can be with or without uncertainty) using stdscore:

stdscore(1.3 ± 0.12, 1)
# => 2.5000000000000004
stdscore(4.7 ± 0.58, 5 ± 0.01)
# => -0.5172413793103445 ± 0.017241379310344827

weightedmean Function¶

Calculate the weighted and arithmetic means of your set of measurements with weightedmean and mean respectively:

weightedmean((3.1±0.32, 3.2±0.38, 3.5±0.61, 3.8±0.25))
# => 3.4665384454054498 ± 0.16812474090663868
mean((3.1±0.32, 3.2±0.38, 3.5±0.61, 3.8±0.25))
# => 3.4000000000000004 ± 0.2063673908348894

TODO¶

• Add pretty printing: optionally print only the relevant significant digits (issue #5)
• Other suggestions welcome :-)

The Measurement Type¶

Measurement is a composite parametric type, whose parameter is the AbstractFloat subtype of the nominal value and the uncertainty of the measurement. Measurement type itself is subtype of AbstractFloat, thus Measurement objects can be used in any function taking AbstractFloat arguments without redefining it, and calculation of uncertainty will be exact.

In detail, this is the definition of the type:

immutable Measurement{T<:AbstractFloat} <: AbstractFloat
    val::T
    err::T
    tag::Float64
    der::Derivatives{Tuple{T, T, Float64}, T}
end

The fields represent:

• val: the nominal value of the measurement
• err: the uncertainty, assumed to be standard deviation
• tag: a unique identifier, it is used to identify a specific measurement in the list of derivatives. This is automatically created with rand. The result of mathematical operation will have this field set to NaN because it is not relevant for non independent measurements.
• der: the list of derivates with respect to the independent variables from which the expression comes. Derivatives is a lightweight dictionary type. The keys are the tuples (val, err, tag) of all independent variables from which the object has been derived, while the corresponding value is the partial derivative of the object with respect to that independent variable.

As already explained in the “Usage” section, every time you use one of the constructors

Measurement(value, uncertainty)
value ± uncertainty

you define a new independent measurement. This happens because these contructors generate a new random and (hopefully) unique tag field, that is used to distinguish between really equal objects and measurements that only by chance share the same nominal value and uncertainty. For these reasons,

x = 24.3 ± 2.7
y = 24.3 ± 2.7

will produce two independent measurements and they will be treated as such when performing mathematical operations. In particular, you can also notice that they are not egal

x === y
# => false

If you instead intend to make y really the same thing as x you have to use assignment:

x = y = 24.3 ± 2.7
x === y
# => true

Thanks to how the Julia language is designed, support for complex measurements, arbitrary precision calculations and array operations came with practically no effort during the development of the package. As explained by Steven G. Johnson, since in Julia a lot of nonlinear functions are internally implemented in terms of elementary operations on the real and imaginary parts it was natural to make the type subtype of Real in order to easily work with complex measurements. In particular, it was then chosen to select the AbstractFloat type because some functions of complex arguments (like sqrt and log) take Complex{AbstractFloat} arguments instead of generic Complex{Real}, and any operation on a Measurement{R} object, with R subtype of Real different from AbstractFloat, would turn it into Measurement{F}, with F subtype of AbstractFloat, anyway.

Correlation¶

One must carefully take care of correlation between two measurements when propagating the uncertainty for an operation. Actually, the term “correlation” may refer to different kind of dependences between two or more quantities, what we mean by this term in Measurements.jl is explained in the “Usage” section of this manual.

Dealing with functional correlation between Measurement objects, when using functions with arity larger than one, is an important feature of this package. This is accomplished by keeping inside each Measurement object the list of its derivatives with respect to the independent variables from which the quantity comes, and this list is updated every time a new mathematical operation is performed. This role is played by the der field. This dictionary is useful in order to trace the contribution of each measurement and propagate the uncertainty in the case of functions with more than one argument.

The use of the list of derivatives has been inspired by Python package uncertainties, but the rest of the implementation of Measurements.jl is completely independent from that of uncertainties package, even though it may happen to be similar.

Uncertainty Propagation¶

For a function \(G(a, b, c, \dots)\) of real arguments with uncertainties \(a = \bar{a} \pm \sigma_{a}\), \(b = \bar{b} \pm \sigma_{b}\), and \(c = \bar{c} \pm \sigma_{c}\), ..., the linear error propagation theory prescribes that uncertainty is propagated as follows:

where the \(\sigma_{ab}\) factors are the covariances defined as

\(E[a]\) is the expected value, or mean, of \(a\). If uncertainties of the quantities \(a\), \(b\), \(c\), ..., are independent and normally distributed, the covariances are null and the above formula for uncertainty propagation simplifies to

In general, calculating the covariances is not an easy task. The trick adopted in Measurements.jl in order to deal with simple functional correlation is to propagate the uncertainty always using really independent variables. Thus, dealing with functional correlation boils down to finding the set of all the independent measurements on which an expression depends. If this set is made up of \(\{x, y, z, \dots\}\), it is possible to calculate the uncertainty of \(G(a, b, c, \dots)\) with

where all covariances due to functional correlation are null. This explains the purpose of keeping the list of derivatives with respect to independent variables in Measurement objects: by looking at the der fields of \(a\), \(b\), \(c\), ..., it is possible to determine the set of independent variables. If other types of correlation (not functional) between \(x\), \(y\), \(z\), ..., are present, they should be treated by calculating the covariances as shown above.

For a function of only one argument, \(G = G(a)\), there is no problem of correlation and the uncertainty propagation formula in the linear approximation simply reads

even if \(a\) is not an independent variable and comes from operations on really independent measurements.

For example, suppose you want to calculate the function \(G = G(a, b)\) of two arguments, and \(a\) and \(b\) are functionally correlated, because they come from some mathematical operations on really independent variables \(x\), \(y\), \(z\), say \(a = a(x, y)\), \(b = b(x, z)\). By using the chain rule, the uncertainty on \(G(a, b)\) is calculated as follows:

What Measurements.jl really does is to calulate the derivatives like \(\partial a/\partial x\) and \(\partial G/\partial x = (\partial G/\partial a)(\partial a/\partial x) + (\partial G/\partial b)(\partial b/\partial x)\), and store them in the der field of \(a\) and \(G\) respectively in order to be able to perform further operations involving these quantities.

Defining Methods for Mathematical Operations¶

Measurements.jl defines new methods for mathematical operations in order to make them accept Measurement arguments. The single most important thing to know about how to define new methods in the package is the Measurements.result. This function, not exported because it is intended to be used only within the package, takes care of propagating the uncertainty as described in the section above. It has two methods: one for functions with arity equal to one, and the other for any other case. This is its syntax:

result(val::Real, der::Real, a::Measurement)

for functions of one argument, and

result(val::Real, der::Tuple{Vararg{Real}},
       a::Tuple{Vararg{Measurement}})

for functions of two or more arguments. The arguments are:

• val: the nominal result of the operation \(G(a, \dots)\);
• der: the partial derivative \(\partial G/\partial a\) of a function \(G = G(a)\) with respect to the argument \(a\) for one-argument functions or the tuple of partial derivatives with respect to each argument in other cases;
• a: the argument(s) of \(G\), in the same order as the corresponding derivatives in der argument.

In the case of functions with arity larger than one, der and a tuples must have the same length.

For example, for a one-argument function like \(\cos\) we have

cos(a::Measurement) = result(cos(a.val), -sin(a.val), a)

Instead, the method for subtraction operation is defined as follows:

-(a::Measurement, b::Measurement) =
    result(a.val - b.val, (1.0, -1.0), (a, b))

Thus, in order to support Measurement argument(s) for a new mathematical operation you have to calculate the result of the operation, the partial derivatives of the functon with respect to all arguments and then pass this information to Measurements.result function.