Examples

These are some basic examples of use of the package:

julia> using Measurements

julia> a = measurement(4.5, 0.1)
4.5 ± 0.1

julia> b = 3.8 ± 0.4
3.8 ± 0.4

julia> 2a + b
12.8 ± 0.4472135954999579

julia> a - 1.2b
-0.05999999999999961 ± 0.49030602688525043

julia> l = measurement(0.936, 1e-3);

julia> T = 1.942 ± 4e-3;

julia> g = 4pi^2*l/T^2
9.797993213510699 ± 0.041697817535336676

julia> c = measurement(4)
4.0 ± 0.0

julia> a*c
18.0 ± 0.4

julia> sind(94 ± 1.2)
0.9975640502598242 ± 0.0014609761696991563

julia> x = 5.48 ± 0.67;

julia> y = 9.36 ± 1.02;

julia> log(2x^2 - 3.4y)
3.3406260917568824 ± 0.5344198747546611

julia> atan2(y, x)
1.0411291003154137 ± 0.07141014208254456

Measurements from Strings

You can construct Measurement objects from strings. Within parentheses there is the uncertainty referred to the corresponding last digits.

julia> measurement("-12.34(56)")
-12.34 ± 0.56

julia> measurement("+1234(56)e-2")
12.34 ± 0.56

julia> measurement("123.4e-1 +- 0.056e1")
12.34 ± 0.56

julia> measurement("(-1.234 ± 0.056)e1")
-12.34 ± 0.56

julia> measurement("1234e-2 +/- 0.56e0")
12.34 ± 0.56

julia> measurement("-1234e-2")
-12.34 ± 0.0

Correlation Between Variables

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

julia> x = 8.4 ± 0.7
8.4 ± 0.7

julia> x - x
0.0 ± 0.0
julia> x/x
1.0 ± 0.0

julia> x*x*x - x^3
0.0 ± 0.0

julia> sin(x)/cos(x) - tan(x)
-2.220446049250313e-16 ± 0.0
# They are equal within numerical accuracy

julia> y = -5.9 ± 0.2

julia> 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:

julia> u = 2x

julia> (x + x) - u
0.0 ± 0.0

julia> u/2x
1.0 ± 0.0

julia> u^3 - 8x^3
0.0 ± 0.0

julia> 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
julia> v = 16.8 ± 1.4

julia> (x + x) - v
0.0 ± 1.979898987322333

julia> v / 2x
1.0 ± 0.11785113019775792
julia> v^3 - 8x^3
0.0 ± 1676.4200705455657

julia> 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.

julia> @uncertain (x -> complex(zeta(x), exp(eta(x)^2)))(2 ± 0.13)
(1.6449340668482273 ± 0.12188127308075564) + (1.9668868646839253 ± 0.042613944993428333)im

julia> @uncertain log(9.4 ± 1.3, 58.8 ± 3.7)
1.8182372640255153 ± 0.11568300475873611

julia> 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

julia> using Cuba

julia> cubaerf(x::Real) =
           2x/sqrt(pi)*cuhre((t, f) -> f[1] = exp(-abs2(t[1]*x)))[1][1]
cubaerf (generic function with 1 method)

julia> @uncertain cubaerf(0.5 ± 0.01)
0.5204998778130466 ± 0.008787825789336267

julia> 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. Remember that if you want to numerically integrate a function that returns a Measurement object you can use QuadGK.jl package, which is written purely in Julia and in addition allows you to set Measurement objects as endpoints, see below.

Tip

Note that the argument of @uncertain macro must be a function call whose arguments are Measurement objects. Thus,

julia> @uncertain zeta(13.4 ± 0.8) + eta(8.51 ± 0.67)

will not work because here the outermost function is +, whose arguments are zeta(13.4 ± 0.8) and eta(8.51 ± 0.67), that however cannot be calculated. You can use the @uncertain macro on each function separately:

julia> @uncertain(zeta(13.4 ± 0.8)) +  @uncertain(eta(8.51 ± 0.67))
1.9974303172187315 ± 0.0012169293212062773

The type of all the arguments provided must be Measurement. If one of the arguments is actually an exact number (so without uncertainty), promote it to Measurement type:

julia> atan2(10, 13.5 ± 0.8)
0.6375487981386927 ± 0.028343666961913202

julia> @uncertain atan2(10 ± 0, 13.5 ± 0.8)
0.6375487981386927 ± 0.028343666962347438

In addition, the function must be differentiable in all its arguments. For example, the polygamma function of order \(m\), polygamma(m, x), is the \(m+1\)-th derivative of the logarithm of gamma function, and is not differentiable in the first argument. Not even the trick of passing an exact measurement would work, because the first argument must be an integer. You can easily work around this limitation by wrapping the function in a single-argument function:

julia> @uncertain (x -> polygamma(0, x))(4.8 ± 0.2)
1.4608477407291167 ± 0.046305812845734776

julia> digamma(4.8 ± 0.2)   # Exact result
1.4608477407291167 ± 0.04630581284451362

Complex Measurements

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

julia> u = complex(32.7 ± 1.1, -3.1 ± 0.2)

julia> v = complex(7.6 ± 0.9, 53.2 ± 3.4)

julia> 2u + v
(73.0 ± 2.3769728648009427) + (47.0 ± 3.4234485537247377)im

julia> sqrt(u * v)
(33.004702573592 ± 1.0831254428098636) + (25.997507418428984 ± 1.1082833691607152)im

You can also verify the Euler’s formula

julia> cis(u)
(6.27781144696534 ± 23.454542573739754) + (21.291738410228678 ± 8.112997844397572)im

julia> cos(u) + sin(u)*im
(6.277811446965339 ± 23.454542573739754) + (21.291738410228678 ± 8.112997844397572)im

Arbitrary Precision Calculations

If you performed an exceptionally good experiment that gave you extremely precise results (that is, with very low relative error), you may want to use arbitrary precision (or multiple precision) calculations, in order not to loose significance of the experimental results. Luckily, Julia natively supports this type of arithmetic and so Measurements.jl does. You only have to create Measurement objects with nominal value and uncertainty of type BigFloat.

Tip

As explained in the Julia documentation, it is better to use the big string literal to initialize an arbitrary precision floating point constant, instead of the BigFloat and big functions. See examples below.

For example, you want to measure a quantity that is the product of two observables \(a\) and \(b\), and the expected value of the product is \(12.00000007\). You measure \(a = 3.00000001 \pm (1\times 10^{-17})\) and \(b = 4.0000000100000001 \pm (1\times 10^{-17})\) and want to compute the standard score of the product with stdscore(). Using the ability of Measurements.jl to perform arbitrary precision calculations you discover that

julia> a = big"3.00000001" ± big"1e-17"

julia> b = big"4.0000000100000001" ± big"1e-17"

julia> stdscore(a * b, big"12.00000007")
7.999999997599999878080000420160000093695993825308195353920411656927305928530607

the measurement significantly differs from the expected value and you make a great discovery. Instead, if you used double precision accuracy, you would have wrongly found that your measurement is consistent with the expected value:

julia> stdscore((3.00000001 ± 1e-17)*(4.0000000100000001 ± 1e-17), 12.00000007)
0.0

and you would have missed an important prize due to the use of an incorrect arithmetic.

Of course, you can perform any mathematical operation supported in Measurements.jl using arbitrary precision arithmetic:

julia> hypot(a, b)
5.000000014000000080399999974880000423919999216953595312794907845334503498479533 ± 1.000000000000000000000000000000000000000000000000000000000000000000000000000009e-17

julia> log(2a) ^ b
1.030668110995484998145373137400169442058573718746529435800255440973153647087416e+01 ± 9.744450581349822034766870718391736028419817951565653507621645979913795265663606e-17

Operations with Arrays and Linear Algebra

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

julia> A = [1.03 ± 0.14, 2.88 ± 0.35, 5.46 ± 0.97]
3-element Array{Measurements.Measurement{Float64},1}:
 1.03±0.14
 2.88±0.35
 5.46±0.97

julia> B = [0.92 ± 0.11, 3.14 ± 0.42, 4.67 ± 0.58]
3-element Array{Measurements.Measurement{Float64},1}:
 0.92±0.11
 3.14±0.42
 4.67±0.58

julia> exp.(sqrt.(B)) .- log.(A)
3-element Array{Measurements.Measurement{Float64},1}:
  2.57996±0.202151
  4.82484±0.707663
  6.98252±1.17829

julia> @. 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

If you originally have separate arrays of values and uncertainties, you can create an array of Measurement objects using measurement or ± with the dot syntax for vectorizing functions:

julia> C = measurement.([174.9, 253.8, 626.3], [12.2, 19.4, 38.5])
3-element Array{Measurements.Measurement{Float64},1}:
 174.9±12.2
 253.8±19.4
 626.3±38.5

julia> sum(C)
1055.0 ± 44.80457565918909

julia> D = [549.4, 672.3, 528.5] .± [7.4, 9.6, 5.2]
3-element Array{Measurements.Measurement{Float64},1}:
 549.4±7.4
 672.3±9.6
 528.5±5.2

julia> mean(D)
583.4 ± 4.396463225012679

Tip

prod and sum (and mean, which relies on sum) functions work out-of-the-box with any iterable of Measurement objects, like arrays or tuples. However, these functions have faster methods (quadratic in the number of elements) when operating on an array of Measurement s than on a tuple (in this case the computational complexity is cubic in the number of elements), so you should use an array if performance is crucial for you, in particular for large collections of measurements.

Some linear algebra functions work out-of-the-box, without defining specific methods for them. For example, you can solve linear systems, do matrix multiplication and dot product between vectors, find inverse, determinant, and trace of a matrix, do LU and QR factorization, etc.

julia> A = [(14 ± 0.1) (23 ± 0.2); (-12 ± 0.3) (24 ± 0.4)]
2×2 Array{Measurements.Measurement{Float64},2}:
  14.0±0.1  23.0±0.2
 -12.0±0.3  24.0±0.4

julia> b = [(7 ± 0.5), (-13 ± 0.6)]
2-element Array{Measurements.Measurement{Float64},1}:
   7.0±0.5
 -13.0±0.6

# Solve the linear system Ax = b
julia> x = A \ b
2-element Array{Measurements.Measurement{Float64},1}:
  0.763072±0.0313571
 -0.160131±0.0177963

# Verify this is the correct solution of the system
julia> A * x ≈ b
true

julia> dot(x, b)
7.423202614379084 ± 0.5981875954418516

julia> det(A)
611.9999999999999 ± 9.51262319236918

julia> trace(A)
38.0 ± 0.4123105625617661

julia> A * inv(A) ≈ eye(A)
true

julia> lufact(A)
Base.LinAlg.LU{Measurements.Measurement{Float64},Array{Measurements.Measurement{Float64},2}} with factors L and U:
Measurements.Measurement{Float64}[1.0±0.0 0.0±0.0; -0.857143±0.0222861 1.0±0.0]
Measurements.Measurement{Float64}[14.0±0.1 23.0±0.2; 0.0±0.0 43.7143±0.672403]

julia> qrfact(A)
Base.LinAlg.QR{Measurements.Measurement{Float64},Array{Measurements.Measurement{Float64},2}}(Measurements.Measurement{Float64}[-18.4391±0.209481 -1.84391±0.522154; -0.369924±0.00730266 33.1904±0.331267],Measurements.Measurement{Float64}[1.75926±0.00836088,0.0±0.0])

Derivative, Gradient and Uncertainty Components

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 a convenient function, Measurements.derivative(), that returns the partial derivative of an expression with respect to independent measurements. Its vectorized version can be used to compute the gradient of an expression with respect to multiple independent measurements.

julia> x = 98.1 ± 12.7
98.1 ± 12.7

julia> y = 105.4 ± 25.6
105.4 ± 25.6

julia> z = 78.3 ± 14.1
78.3 ± 14.1

julia> Measurements.derivative(2x - 4y, x)
2.0

julia> Measurements.derivative(2x - 4y, y)
-4.0

julia> Measurements.derivative.(log1p(x) + y^2 - cos(x/y), [x, y, z])
3-element Array{Float64,1}:
   0.0177005
 210.793
   0.0       # The expression does not depend on z

Tip

The vectorized version of Measurements.derivative() is useful in order to discover which variable contributes most to the total uncertainty of a given expression, if you want to minimize it. This can be calculated as the Hadamard (element-wise) product between the gradient of the expression with respect to the set of variables and the vector of uncertainties of the same variables in the same order. For example:

julia> w = y^(3//4)*log(y) + 3x - cos(y/x)
447.0410543780643 ± 52.41813324207829

julia> abs.(Measurements.derivative.(w, [x, y]) .* Measurements.uncertainty.([x, y]))
2-element Array{Float64,1}:
 37.9777
 36.1297

In this case, the x variable contributes most to the uncertainty of w. In addition, note that the Euclidean norm of the Hadamard product above is exactly the total uncertainty of the expression:

julia> vecnorm(Measurements.derivative.(w, [x, y]) .* Measurements.uncertainty.([x, y]))
52.41813324207829

The Measurements.uncertainty_components() function simplifies calculation of all uncertainty components of a derived quantity:

julia> Measurements.uncertainty_components(w)
Dict{Tuple{Float64,Float64,Float64},Float64} with 2 entries:
  (98.1, 12.7, 0.303638)  => 37.9777
  (105.4, 25.6, 0.465695) => 36.1297

julia> vecnorm(collect(values(Measurements.uncertainty_components(w))))
52.41813324207829

stdscore Function

You can get the distance in number of standard deviations between a measurement and its expected value (not a Measurement) using stdscore():

julia> stdscore(1.3 ± 0.12, 1)
2.5000000000000004

You can use the same function also to test the consistency of two measurements by computing the standard score between their difference and zero. This is what stdscore() does when both arguments are Measurement objects:

julia> stdscore((4.7 ± 0.58) - (5 ± 0.01), 0)
-0.5171645175253433

julia> stdscore(4.7 ± 0.58, 5 ± 0.01)
-0.5171645175253433

weightedmean Function

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

julia> weightedmean((3.1±0.32, 3.2±0.38, 3.5±0.61, 3.8±0.25))
3.4665384454054498 ± 0.16812474090663868

julia> mean((3.1±0.32, 3.2±0.38, 3.5±0.61, 3.8±0.25))
3.4000000000000004 ± 0.2063673908348894

Measurements.value and Measurements.uncertainty Functions

Use Measurements.value() and Measurements.uncertainty() to get the values and uncertainties of measurements. They work with real and complex measurements, scalars or arrays:

julia> Measurements.value(94.5 ± 1.6)
94.5

julia> Measurements.uncertainty(94.5 ± 1.6)
1.6

julia> Measurements.value.([complex(87.3 ± 2.9, 64.3 ± 3.0), complex(55.1 ± 2.8, -19.1 ± 4.6)])
2-element Array{Complex{Float64},1}:
 87.3+64.3im
 55.1-19.1im

julia> Measurements.uncertainty.([complex(87.3 ± 2.9, 64.3 ± 3.0), complex(55.1 ± 2.8, -19.1 ± 4.6)])
2-element Array{Complex{Float64},1}:
 2.9+3.0im
 2.8+4.6im

Interplay with Third-Party Packages

Measurements.jl works out-of-the-box with any function taking arguments no more specific than AbstractFloat. This makes this library particularly suitable for cooperating with well-designed third-party packages in order to perform complicated calculations always accurately taking care of uncertainties and their correlations, with no effort for the developers nor users.

The following sections present a sample of packages that are known to work with Measurements.jl, but many others will interplay with this package as well as them.

Numerical Integration with QuadGK.jl

The powerful integration routine quadgk from QuadGK.jl package is smart enough to support out-of-the-box integrand functions that return arbitrary types, including Measurement:

julia> QuadGK.quadgk(x -> exp(x / (4.73 ± 0.01)), 1, 7)
(14.933307243306032 ± 0.009999988180463411, 0.0 ± 0.010017961523508253)

Measurements.jl pushes the capabilities of quadgk further by supporting also Measurement objects as endpoints:

julia> QuadGK.quadgk(cos, 1.19 ± 0.02, 8.37 ± 0.05)
(-0.05857827689796702 ± 0.02576650561689427, 2.547162480937004e-11)

Compare this with the expected result:

julia> sin(8.37 ± 0.05) - sin(1.19 ± 0.02)
-0.058578276897966686 ± 0.02576650561689427

Also with quadgk correlation is properly taken into account:

julia> a = 6.42 ± 0.03
6.42 ± 0.03

julia> QuadGK.quadgk(sin, -a, a)
(2.484178227707412e-17 ± 0.0, 0.0)

If instead the two endpoints have, by chance, the same nominal value and uncertainty but are not correlated:

julia> QuadGK.quadgk(sin, -6.42 ± 0.03, 6.42 ± 0.03)
(2.484178227707412e-17 ± 0.005786464233000303, 0.0)

Numerical and Automatic Differentiation

With Calculus.jl package it is possible to perform numerical differentiation using finite differencing. You can pass in to the Calculus.derivative function both functions returning Measurement objects and a Measurement as the point in which to calculate the derivative.

julia> using Measurements, Calculus

julia> a = -45.7 ± 1.6
-45.7 ± 1.6

julia> b = 36.5 ± 6.0
36.5 ± 6.0

julia> Calculus.derivative(exp, a) ≈ exp(a)
true

julia> Calculus.derivative(cos, b) ≈ -sin(b)
true

julia> Calculus.derivative(t -> log(-t * b)^2, a) ≈ 2log(-a * b)/a
true

Other packages provide automatic differentiation methods. Here is an example with AutoGrad.jl, just one of the packages available:

julia> using AutoGrad

julia> grad(exp)(a) ≈ exp(a)
true

julia> grad(cos)(b) ≈ -sin(b)
true

julia> grad(t -> log(-t * b)^2)(a) ≈ 2log(-a * b)/a
true

However remember that you can always use Measurements.derivative() to compute the value (without uncertainty) of the derivative of a Measurement object.

Use with SIUnits.jl and Unitful.jl

You can use Measurements.jl in combination with a third-party package in order to perform calculations involving physical measurements, i.e. numbers with uncertainty and physical unit. The details depend on the specific package adopted. Such packages are, for instance, SIUnits.jl and Unitful.jl. You only have to use the Measurement object as the value of the SIQuantity object (for SIUnits.jl) or of the Quantity object (for Unitful.jl). Here are a few examples.

julia> using Measurements, SIUnits, SIUnits.ShortUnits

julia> hypot((3 ± 1)*m, (4 ± 2)*m) # Pythagorean theorem
5.0 ± 1.7088007490635064 m

julia> (50 ± 1)Ω * (13 ± 2.4)*1e-2*A # Ohm's Law
6.5 ± 1.20702112657567 kg m²s⁻³A⁻¹

julia> 2pi*sqrt((5.4 ± 0.3)*m / ((9.81 ± 0.01)*m/s^2)) # Pendulum's  period
4.661677707464357 ± 0.1295128435999655 s


julia> using Measurements, Unitful

julia> hypot((3 ± 1)*u"m", (4 ± 2)*u"m") # Pythagorean theorem
5.0 ± 1.7088007490635064 m

julia> (50 ± 1)*u"Ω" * (13 ± 2.4)*1e-2*u"A" # Ohm's Law
6.5 ± 1.20702112657567 A Ω

julia> 2pi*sqrt((5.4 ± 0.3)*u"m" / ((9.81 ± 0.01)*u"m/s^2")) # Pendulum's period
4.661677707464357 ± 0.12951284359996548 s