/
05-project-errorprop.jmd
139 lines (97 loc) · 6.5 KB
/
05-project-errorprop.jmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
---
title : "Project: automatic error propagation"
author : Michiel Stock, Bram De Jaegher, Daan Van Hauwermeiren
date : 18 December 2019
---
# What is error propagation?
During my first year of college, I really liked the physics lectures because everything was exact. In my first year of physics, I really disliked the physics labs, because everything was messy. First-year physics labs are generally disliked partially because of the uninspiring topics (measuring resistance in a wire! determining a heat transfer coefficient! measuring chirality of sugar!) and partly because this was the only lab in which 'exact' measurements had to be performed. These labs introduced the complex and unexplained rules for measurement error propagation. Without any prior statistics or probability courses, it was lost on me where these strange rules originated from, as sharp contrast to nearly everything else in physics.
A year later, after plowing through a basic probability course, error propagation is less mysterious. Almost all these rules to account for the error can be derived from a simple principle:
> Given two **independent** random variables $$X$$ and $$Y$$, the variance of a linear combination $$\text{Var}[aX+bY]=a^2\text{Var}[X] + b^2\text{Var}[Y]$$.
Measurement errors are usually given by a standard deviation, the square root of the variance. Given this principle, error propagation is merely bookkeeping of the standard error on the measurement for various computations. This [table on error propagation](https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Example_formulae) might be useful.
For nonlinear functions, we can compute an approximate uncertainty propagation using a first-order Taylor approximation. We have, for any function $f(x)$:
$$
f(x\pm \sigma) \approx f(x) \pm |f'(x)|\sigma\,.
$$
For example, for squaring a function, we have
$$
(x\pm\sigma)^2 = x^2 \pm 2|x|\sigma\,.
$$
Note that this is consistent with the above rules for multiplication. Let us implement the general formula for raising a measurement to the power $p$.
We can implement this for all the standard mathematical functions one by one. However, Julia provides use with two efficient tools to do this in one swoop: automatic differentiation and metaprogramming. We just loop a list of functions of interest and automatically generate the correct approximate rule.
Instead of processing the measurements and the standard error separately, suppose we could just make a new type of number which contains both the observed value and its uncertainty. And suppose we could just compute something with these numbers, plugging them into our formulas where the error is automatically accounted for using the standard error propagation rules. In Julia, it is dead simple to construct such new numbers and just overload existing functions such that they compute in the correct way.
# Goal of this project
We will implement a binary operator `±`, which can be used to add the standard error to a measurement.
The result is a new type `Measurement` containing both the value and the standard error. Standard functions will be overloaded to process this structure correctly.
# Assignments
1. Make an `Measurement` structure with two fields: the measured value `x` and its standard deviation `σ` (`\sigma<TAB>`). Make sure that both `x` and `σ` are of the same type and a subtype of `Real`.
2. Make a constructor for `Measurement`, such that an error is returned when a negative standard error is given.
3. Make two functions `val` and `err`, which respecitivly return the value and the standard error of a measurement.
4. Make a binary operator `±` (`\pm<TAB>`), such that `x ± σ` returns a new instance of the type `Measurement`. Try it on a value 4.0 with a standard error of 1.2.
5. Overload all functions of scalar multiplication, adding and substracting measurements and adding a constant such that they correctly process measurements. Try some examples.
6. Overload the power function `^` such that one can raise `Measurement`s to a power. Note the special case `(x ± σ)^0` = `one(x) ± zero(σ)`.
7. Run example to generate error functions for a whole range of functions using metaprogramming.
8. Solve the small exercise using data.
**Assignment 1 and 2**
```julia; eval=false
struct Measurement
...
function Measurement(...)
...
end
end
```
**Assignment 3**
```julia; eval=false
val(m::Measurement) = ...
err(m::Measurement) = ...
```
**Assignment 4**
Don't forget to add type annotations!
```julia; eval=false
±(x, σ) = ...
```
**Assignment 5**
```julia; eval=false
# scalar multiplication
Base.:*(a::Real, m::Measurement) = ...
Base.:/(m::Measurement, a::Real) = ...
# adding and substracting measurements
Base.:+(m1::Measurement, m2::Measurement) = ...
Base.:-(m1::Measurement, m2::Measurement) = ...
Base.:-(m::Measurement) = ...
# adding a constant
Base.:+(m::Measurement, a::Real) = ...
Base.:+(a::Real, m::Measurement) = ...
# multiplying two measurments
Base.:*(m1::Measurement, m2::Measurement) = ...
```
**Assignment 6**
```julia; eval=false
Base.:^(m::Measurement, p::Integer) = ...
```
**Assignment 7**
```julia; eval=false
using ForwardDiff
for f in [:sin, :cos, :tan, :exp, :log, :log2, :log10, :sqrt, :inv]
eval(quote
# this is a line of code generated using string interpolation
Base.$f(m::Measurement) = $f(m.x) ± abs(ForwardDiff.derivative($f, m.x) * m.σ)
end)
end
```
**Assignment 8**
Let's apply this to in a somewhat realistic setting. Many methods in analytical chemistry are based on the law of [Beer-Lambert](https://en.wikipedia.org/wiki/Beer%E2%80%93Lambert_law). This law relates the absorption of a ray of light passing through a cuvet with the concentration of a solution. For a given reference intensity $I_0$ at concentration of 0 and a lower intensity $I$ when it passes through a solution of concentration $c$, we have
$$
\log \left(\frac{I_0}{I}\right) = \varepsilon c l\,,
$$
with $\varepsilon$ the molar extinction coefficient and $l$ the thickness of the cuvet.
Suppose we want to determine the extinction coefficient for some substance, using a cuvet of thickness $0.02\pm 0.001$ and a reference solution of a concentration of $0.73\pm 0.02$ M.
We perform some intensity measurements, with associated measurement errors.
| $I0$ | $\sigma_{I0}$ | $I$ | $\sigma_{I}$ |
|------|---------------|------|---------------|
| 0.8 | 0.14 | 0.2 | 0.12 |
| 1.1 | 0.11 | 0.3 | 0.07 |
| 1.2 | 0.08 | 0.4 | 0.101 |
**Estimate the molecular extinction coefficient.**
```julia
```