Transcript Document

Some Advanced Topics using Propagation of Errors and Least Squares Fitting
More on Least Squares Fit (LSQF)
In Lec 5, we discussed how we can fit our data points to a linear function (straight line) and
get the "best" estimate of the slope and intercept.
However, we did not discuss two important issues:
I) How to estimate the uncertainties on our slope and intercept obtained from a LSQF?
II) How to apply the LSQF when we have a non-linear function?
Estimation of Errors on parameters determined from a LSQF
Assume we have data points that lie on a straight line:
y = a + bx
Assume we have n measurements of the x’s and y's.
For simplicity, assume that each y measurement has the same error, .
Assume that x is known much more accurately than y.
 ignore any uncertainty associated with x.
Previously we showed that the solution for the intercept a and slope b is:
n
a
R. Kass/Sp07
n
n
n
2
 yi  xi   xi yi  xi
i1 i1
i1
i1
n
n
n  xi2  (  xi )2
i1
i1
and b 
n
n
n
i1
n
i1
n
i1
n  xi yi   xi  yi
n  xi2  (  xi )2
i1
i1
P416/Lecture 7
1
Since a and b are functions of the measurements (yi's) we can use the Propagation of Errors technique to
estimate a and b.
 Q2
2
2
Q Q 
2 Q 
2 Q 
  x     y    2 xy   
See lecture 4
x 
(general formula for 2 variables)
x y 
y 
Assume that: a) Each measurement is independent of each other (xy=0).
b) We can neglect any error () associated with x.
c) Each y measurement has the same , i.e. i=.

 Q2
 Q 
  x2  
2
 Q 
  y2  
2
 Q2
 x 
a 
2
2
assumption a)
 y 
assumption b)
 y 
2  a
  yi 
i 1
 yi
n
n
a

yi
 Q 
  y2  
2
n  a 

   2   
i 1 yi 

2
n
n
2
y x  x y x
 i 1 i j 1 j i 1 i i j 1 j
n
n
yi
n  xi2  (  xi ) 2
i 1
i 1
assumption c)
n
n

n
2
 x j  xi  x j
j 1
j 1
n
n
n  xi2  (  xi ) 2
i 1
i 1
plugging our formula for a
2
 n x2  x n x 
 ( n x 2 )2  x 2 ( n x )2  2 x n x n x 2 



  j
j
i
j 
i  j
i  j  j 
n
n
j 1
j 1
j 1 j 1
   2   j 1

 a2   2   j n1
n
n
n
2
2

i 1
i 1
(n  xi2  (  xi ) 2 ) 2
 n  xi  (  xi ) 


i 1
i 1
i 1
 i 1



R. Kass/Sp07
P416/Lecture 7
2
 a2   2
n
n
n
n
n
j1
i1
j1
j1
j1
n(  x 2j )2   xi2 (  x j )2  2(  x j )2  x 2j
n
n
2
(n  xi  (  xi )2 )2
i1
i1
n
n
2
n  x j  (  x j )2
n
j1
2
   x 2j j1
n
n
j1 (n x 2  ( x )2 )2
 i
 i
i1
i1
n
2
xj
 a2   2 n j1 n
variance
2
2
n  xi  (  xi )
i1
i1
2
n
n
n
j1
i1
n
j1
n(  x 2j )2   xi2 (  x j )2
n
(n  xi2  (  xi )2 )2
i1
i1
in the intercept
We can find the variance in the slope (b) using exactly the same procedure:

n

2
2
2
 nxi   x j 
n
n
n




b
b
j1

 b2    2y i     2      2   n
n
yi 
i1
i1yi 
i1n x 2  ( x )2 
 i 
  i
i1
 i1

n
2
R. Kass/Sp07
2
n
n
n
n
2
2
 x j  n(  x j )  2n  xi  x j
j1
j1
i1 j1
n
n
(n  xi2  (  xi )2 )2
i1
i1
P416/Lecture 7
3
n 2
2
b 
n
n
n
n  xi2
i1
n
 (  xi )
variance in the slope
2
i1
If we don't know the true value of , we can estimate the variance using the spread between the
measurements (yi’s) and the fitted values of y:

1 n
1 n
2
fit 2
 
 (yi  yi ) 
 (yi  a  bxi )2
n  2 i1
n  2 i1
n  2 = number of degree of freedom
= number of data points  number of parameters (a,b) extracted from the data
If each yi measurement has a different error i:

1 n xi2
2
a   2
D i1 i
 b2 
weighted slope
and intercept
1 n 1

D i1 i2
n
D 
1
n

xi2
2
2
i1 i i1 i
n
 (
xi
2
)
2
i1 i
The above expressions simplify to the “equal variance” case.
Don't forget to keep track of the “n’s” when factoring out equal ’s. For example:

n

1
2
i1 i
R. Kass/Sp07

n
2
not
1
2
P416/Lecture 7
4
Example: Fit to straight line
x
1.0
2.0
3.0
4.0
y
2.2
2.9
4.3
5.2

0.2
0.4
0.3
0.1
(same data as Lec 6)
Previous fit to: f(x,b)=1+bx, b=1.05
New fit to: f(x,b)=a+bx,
a =1.158±0.253 b=1.011± 0.072
4
a and b calculated using weighted fit
n
n
(Taylor, P 198)
xi2 n yi
xi
a
 
i 1

2
i
i 1

xi yi

2
i i 1
2
i
n
yi
D
n
b
2
i i 1
y
n
1
n
xi yi
  
2
i i 1
i 1
2
i
n

i 1
xi


2
i i 1
D
n
D
i 1
1

n
xi2

2
i i 1
2
i
n
 (
i 1
xi

2
i
1 n xi2
a   2
D i 1  i
2
R. Kass/Sp07
2
1
2
x
3
4
)2
Errors on a and b calculated using:
 b2 
2
i
1 n 1

D i 1  i2
χ2=0.65 for 2 degrees of freedom
P(χ2>0.65 for 2dof) = 72%
P416/Lecture 7
5
l
LSQF with non-linear functions:
u
For our purposes, a non-linear function is a function where one or more of the parameters that
we are trying to determine (e.g. a, b from the straight line fit) is raised to a power other than 1.
n Example: functions that are non-linear in the parameter :
y  A  x /
y  A  x 2
y  Aex / 
However, these functions are linear in the parameters A.
u The problem with most non-linear functions is that we cannot write down a solution for
the parameters in a closed form using, for example, the techniques of linear algebra (i.e. matrices).
 n Usually non-linear problems are solved numerically using a computer.
n Sometimes by a change of variable(s) we can turn a non-linear problem into a linear one.
Example: take the natural log of both sides of the above exponential equation:
lny=lnA-x/=C+Dx
u
= -1/D
Now a linear problem in the parameters C and D!
In fact it’s just a straight line!
To measure the lifetime  (Lab 7) we first fit for D and then transform D into .
Example: Decay of a radioactive substance. Fit the following data to find N0 and :
N (t )  N0et /
n
n
n
N represents the amount of the substance present at time t.
N0 is the amount of the substance at the beginning of the experiment (t = 0).
 is the lifetime of the substance.
R. Kass/Sp07
P416/Lecture 7
6
Lifetime Data
i
1
2
3
4
5
6
7
8
9
10
ti
0
15
30
45
60
75
90
105
120
135
Ni
106
80
98
75
74
73
49
38
37
22
yi = ln Ni
4.663
4.382
4.585
4.317
4.304
4.290
3.892
3.638
3.611
3.091
D
n
n
n
i 1
n
i 1
n
i 1
n  xi yi   yi  xi
n  xi2  (  xi ) 2
i 1

10  2560.41  40.773  675
 0.01033
2
10  64125  (675)
i 1
  1 / D  96.80 seconds
n
The intercept is given by: C = 4.77 = ln A or A = 117.9
5
y = 4.7746 + -0.010331x R= 0.93518
Y(x)
ln[N(t)]
4.5
4
3.5
3
-20
R. Kass/Sp07
0
20
40
60
time, t
x
80
100
120
140
P416/Lecture 7
7
Example: Find the values A and  taking into account the uncertainties in the data points.
n The uncertainty in the number of radioactive decays is governed by Poisson statistics.
n The number of counts Ni in a bin is assumed to be the average (m) of a Poisson distribution:
m = Ni = Variance
n The variance of yi(= ln Ni) can be calculated using propagation of errors:
2
2
2
 2y   N2 y/N   (N) ln N /N   (N)1/ N  1/ N
n The slope and intercept from a straight line fit that includes uncertainties in the data points:
n y n x2
n x y n x
n 1 n x y
n x n y
i
i 
i i
i
i i 
 2 2  2  2
 2  2  i2  i2








Taylor P198 and
a  i1 i i1 i 2 i1 i i1 i and b  i1 i i1 i 2 i1 i i1 i

n 1 n x
n x
n 1 n x
n x
problem 8.9
2
i (
i )2
i (
 2 2  2
 2  2  i2 )
u
i1 i i1 i
i1 i
i1 i i1 i
i1 i
If all the 's are the same then the above expressions are identical to the unweighted case.
a=4.725 and b = -0.00903
 = -1/b = -1/0.00903 = 110.7 sec

n
To calculate the error on the lifetime (), we first must calculate the error on b:
n 1
 2
652
i1 i
6
 b2 


1.0110
n 1 n x2
n x
652  2684700  (33240)2
i (
i )2
 2 2  2
i1 i i1 i
i1 i
2
2
1.005 10 3

 12.3
3 2
(9.03 10 )
 
    b  / b       b 1/ b
2
2
The experimentally determined lifetime is:  = 110.7 ± 12.3 sec.

R. Kass/Sp07
P416/Lecture 7
8
Least Squares Fitting-Exponential Example
We can calculate the c2 to see how “good” the data fits an exponential decay
distribution:
For this problem: lnA=4.725  A=112.73 and  = 110.7 sec
c 2  i 1
10
( N i  Aeti /  ) 2
 i2
( N i  Aeti /  ) 2
 i 1
 15.6
t i / 
Ae
10
Poisson approximation
Mathematica Calculation:
Do[{csq=csq+(cnt[i]-a*Exp[-x[i]/tau])^2/(a*Exp[-x[i]/tau])},{i,1,10}]
Print["The chi sq per dof is ",csq/8]
xvt=1-CDF[ChiSquareDistribution[8],csq];
The chi sq per dof is 1.96
Print["The chi sq prob. is ",100*xvt,"%"]
The chi sq prob. is 4.9 %
This is not such a good fit since the probability is only ~4.9%.
R. Kass/Sp07
P416/Lecture 7
9
The Error on the Mean
Error on the mean (review from Lecture 4)
Question: If we have a set of measurements of the same quantity:
x1  1 x2   2 ... xn   n
What's the best way to combine these measurements?
How to calculate the variance once we combine the measurements?
Assuming Gaussian
statistics, the Maximum Likelihood Method says combine the measurements as:
n

 xi /  i2
x  i1n
weighted average
2
1/  i
i1
If all the variances (21= 22= 23..) are the same:
1 n
x   xi
unweighted average
n i1

The variance of the weighted average can be calculated using propagation of errors:
2
1 /  i2

n



x n
And the derivatives are:
 x2   
x   i2

x

i
i 1  xi

1 /  i2

i 1  xi
n
 x2   
 x2 
2

1/ 
1
2
x   i2  


i
2
2
n
i 1 
n

2
2
 1 /  i 
1 /  i 
 i 1

 i 1

n
4
i
i 1
n
1 / 
i 1
2
i

1
n
1

i 1
2
i
1
n
1

i 1
2
i
R. Kass/Sp07
x=error in the weighted mean
P416/Lecture 7
10
If all the variances are the same:
 x2

1
n
2
1 /  i
i 1

2
n
Lecture 4
The error in the mean (x) gets smaller as the number of measurements (n) increases.
Don't confuse the error in the mean (x) with the standard deviation of the distribution ()!
If we make more measurements:
The standard deviation () of the (gaussian) distribution remains the same,
BUT the error in the mean (x) decreases
Example: Two experiments measure the mass of the proton:
Experiment 1 measures mp=950±25 MeV
Experiment 2 measures mp=936±5 MeV
Using just the average of the two experiments we find:
mp=(950+936)/2=943MeV
Using the weighted average of the two experiments we find:
mp=(950/252+936/52)/(1/252 +1/52)=936.5MeV
and the variance:
2=1/(1/252 +1/52)=24 MeV2
=4.9 MeV
Since experiment 2 was more precise than experiment 1 we would expect the final result to be
closer to experiment 2’s value.
mp=(936.5±4.9) MeV
R. Kass/Sp07
P416/Lecture 7
11