Examples are taken from http://ipl.physics.harvard.edu/wp-uploads/2013/03/PS3_Error_Propagation_sp13.pdf and used on MCMC to show how the answers work
Example: suppose you measure the height H of a door and get 2.00 ± 0.03 m. This means that H = 2.00 m and δH = 0.03 m. The door has a knob which is a height h = 0.88 ± 0.04 m from the bottom of the door. Then the distance from the doorknob to the top of the door is Q = H − h = 1.12 m. What is the uncertainty in Q?
Q = 1.12 ± 0.05 m
In [1]:
import numpy as np
import pymc3 as pm
import seaborn as sns
sns.set(font_scale=1.5)
%matplotlib inline
This version uses prior distributions to do all the work. H and h are both informative priors that then drive the solution to the right answer.
In [2]:
with pm.Model() as model:
H = pm.Normal('H', 2.00, sigma=0.03)
h = pm.Normal('h', 0.88, sigma=0.04)
Q = pm.Deterministic('Q', H-h)
trace = pm.sample(10000)
In [3]:
pm.summary(trace).round(3)
Out[3]:
In [4]:
pm.traceplot(trace, combined=False)
print("MCMC gives {:.2f} +/- {:.2f}, analytic gives {} +/- {}".format(trace['Q'].mean(),
trace['Q'].std(), 1.12, 0.05))
Example: a bird flies a distance d = 120 ± 3 m during a time t = 20.0 ± 1.2 s. The average speed of the bird is v = d/t = 6 m/s. What is the uncertainty of v?
0.39 m/s.
In [5]:
with pm.Model() as model:
d = pm.Normal('d', 123, tau=(3)**-2)
t = pm.Normal('t', 20.0, tau=(1.2)**-2)
v = pm.Deterministic('v', d/t)
trace = pm.sample(40000, chains=4)
In [14]:
pm.summary(trace).round(3)
Out[14]:
In [15]:
pm.traceplot(trace, combined=False, lines=[('d', {}, 123), ('t', {}, 20), ('v', {}, 6)])
print("MCMC gives {0:.2f}, analytic gives {1}".format(trace['v'].std(), 0.39))
Example: the period of an oscillation is measured to be T = 0.20 ± 0.01 s. Thus the frequency is f = 1/T = 5 Hz. What is the uncertainty in f? Answer: the percent uncertainty in T was 0.01/0.20 = 5%. Thus the percent uncertainty in f is also 5%, which means that δf = 0.25 Hz. So f = 5.0 ± 0.3 Hz (after rounding).
f = 5.0 ± 0.3 Hz
In [18]:
with pm.Model() as model:
T = pm.Normal('T', 0.20, tau=(0.01)**-2)
pm.Deterministic('1/T', 1/T)
trace = pm.sample(10000, tune=1000)
In [19]:
pm.traceplot(trace, combined=False)
Out[19]:
In [21]:
pm.summary(trace).round(3)
Out[21]:
In [22]:
print("MCMC gives {0:.1f} +/- {1:.1f}, analytic gives {2} +/- {3}".format(np.mean(trace['1/T']),
np.std(trace['1/T']),
5.0, 0.3))
Example: a ball is tossed straight up into the air with initial speed v0 = 4.0 ± 0.2 m/s. After a time t = 0.60±0.06 s, the height of the ball is y = v0t− 1 2 gt2 = 0.636 m. What is the uncertainty of y? Assume g = 9.80 m/s2 (no uncertainty in g).
Thus y would be properly reported as 0.6 ± 0.4 m.
In [24]:
with pm.Model() as model:
g = 9.80
t = pm.Normal('t', 0.60, tau=(0.06)**-2)
v0 = pm.Normal('v0', 4.0, tau=(0.2)**-2)
h = pm.Deterministic('h', v0*t - 0.5*g*t**2)
trace = pm.sample(10000)
In [25]:
pm.traceplot(trace, combined=False)
Out[25]:
In [26]:
pm.summary(trace).round(3)
Out[26]:
In [27]:
print("MCMC gives {0:.1f} +/- {1:.1f}, analytic gives {2} +/- {3}".format(np.mean(trace['h']),
np.std(trace['h']),
0.6, 0.4))
In [34]:
with pm.Model() as model:
A = pm.Normal('A', 3.6, tau=(0.2)**-2)
B = pm.Normal('B', 3.3, tau=(0.3)**-2)
D = pm.Deterministic('D', A-B)
trace = pm.sample(1000, chains=6)
In [35]:
pm.summary(trace).round(3)
Out[35]:
In [36]:
pm.traceplot(trace, combined=False);
In [37]:
print("MCMC gives {0:.1f} +/- {1:.1f}, analytic gives {2} +/- {3}".format(np.mean(trace['D']),
np.std(trace['D']),
0.3, 0.4))
In [ ]: