{VERSION 4 0 "IBM INTEL NT" "4.0" }
{USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0
1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0
0 0 1 }{CSTYLE "2D Input" 2 19 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 0 14 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
257 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 263 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 268 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 270 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 273 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 275 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
277 "" 1 14 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 278 "" 1 14 0
0 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 279 "" 1 14 0 0 0 0 0 0 0 0
0 0 0 0 0 1 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0
0 0 0 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }{PSTYLE "Heading
1" 0 3 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }1 0 0
0 8 4 0 0 0 0 0 0 -1 0 }{PSTYLE "R3 Font 0" -1 256 1 {CSTYLE "" -1 -1
"Helvetica" 1 18 0 0 0 0 2 1 2 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0
0 0 -1 0 }{PSTYLE "R3 Font 2" -1 257 1 {CSTYLE "" -1 -1 "Courier" 1
11 0 0 0 0 2 1 2 0 0 0 0 0 0 1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }
{PSTYLE "" 0 258 1 {CSTYLE "" -1 -1 "" 1 16 0 0 0 0 0 0 0 0 0 0 0 0 0
1 }0 0 0 -1 -1 -1 0 0 0 0 0 0 -1 0 }}
{SECT 0 {EXCHG {PARA 258 "" 0 "" {TEXT 256 13 "Data analysis" }}{PARA
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 424 "In this worksh
eet we provide some means to analyze experimental data. It represents
an exercise in making use of drawing tools (to plot error bars), in \+
reading in data from an ASCII file (e.g., from a computer-interfaced \+
experiment, spreadsheet/graphing program, etc.), and finally in makin
g use of the statistical analysis package. Ther is also a hand-writte
n least-squares fitting procedure is attached at the end. " }}{PARA
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 198 "First we make \+
use of the data interface. The file draw1 has to exist in the current
directory. We assume four columns of data, containing x_i, y_i, and \+
low and high values for y_i (uncertainties)." }}{PARA 0 "" 0 "" {TEXT
-1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 164 "We follow the conventi
on on PC's for the directory structure. Note, that instead of the forw
ard-slash (borrowed from unix), one could use the double-backward slas
h." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 84 "### WARNING: persistent store
makes one-argument readlib obsolete\nreadlib(readdata):" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "#data:=readdata(`C:\\\\ComPhys\\\\D
ataAnalysis\\\\draw1`,float,4): # this works as well" }}}{EXCHG {PARA
0 "> " 0 "" {MPLTEXT 1 0 56 "data:=readdata(`C:/ComPhys/DataAnalysis/d
raw1`,float,4):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 53 "After reading \+
the data they are available as follows:" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 16 "data[1],data[2];" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 26 "dline:=map(evalf,data[1]);" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 126 "Unfortunately the Maple development team is not listen
ing to user's complaints about the meagre selection of plotting symbol
s." }}{PARA 0 "" 0 "" {TEXT -1 220 "Experimental physicists with data \+
plus uncertainties are the last thing on their minds... So we write ou
r own graphiscs primitives here (they are not very fast when it comes \+
to plotting many points, but one can survive)." }}}{SECT 0 {PARA 3 ""
0 "" {TEXT 257 55 "Graphics primitives to draw data points with error \+
bars" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 86 "To plot a decent-sized ope
n symbol we define our own circle with the centre location [" }{TEXT
259 1 "x" }{TEXT -1 2 ", " }{TEXT 258 1 "y" }{TEXT -1 60 "] passed as \+
the first argument, and the radius as the third." }}{PARA 0 "" 0 ""
{TEXT -1 48 "A fourth argument can be used to set the colour." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "Circle:=proc() local col;" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "if nargs<4 then col:=black else co
l:=args[4] fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "plot([args[3]*cos
(t)+args[1],args[3]*sin(t)+args[2],t=0..2*Pi],scaling" }}{PARA 0 "> "
0 "" {MPLTEXT 1 0 42 "=CONSTRAINED,colour=col,thickness=2); end:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "P1:=Circle(2,3,0.3):" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 26 "A line connecting points [" }{TEXT 260 1
"x" }{TEXT -1 3 "1, " }{TEXT 261 1 "y" }{TEXT -1 5 "1], [" }{TEXT 263
1 "x" }{TEXT -1 3 "2, " }{TEXT 262 1 "y" }{TEXT -1 8 "2] with " }
{TEXT 264 1 "x" }{TEXT -1 3 "1, " }{TEXT 265 1 "y" }{TEXT -1 2 "1," }
{TEXT 266 2 " x" }{TEXT -1 3 "2, " }{TEXT 267 1 "y" }{TEXT -1 99 "2 pa
ssed as the first four arguments to the procedure. The fifth argument \+
allows to set the colour." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
23 "Line:=proc() local col;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "if n
args<5 then col:=black else col:=args[5] fi;" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 70 "plot([[args[1],args[2]],[args[3],args[4]]],scaling=CO
NSTRAINED,colour=" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "col); end:" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "P2:=Line(-1,-1,2,3,red):" }
}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "display(\{P1,P2\},axes=box
ed);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 128 "Here comes the procedure
to plot a data point with error bar. Five arguments specify the point
, and the sixth is for the colour." }}{PARA 0 "" 0 "" {TEXT -1 80 "The
first two arguments provide the location of the point, the next two t
he low-" }{TEXT 269 1 "y" }{TEXT -1 10 " and high-" }{TEXT 268 1 "y" }
{TEXT -1 112 " values, and the fifth argument specifies the width of t
he error bar to either side, and the size of the symbol." }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "Expoint:=proc() local size,col; siz
e:=2*args[5];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "if nargs<6 then co
l:=black else col:=args[6] fi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "C
ircle(args[1],args[2],0.5*args[5],col),Circle(args[1],args[2],args[5]
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 70 ",col),Line(args[1],args[3],args
[1],args[4],col),Line(args[1]-size,args" }}{PARA 0 "> " 0 "" {MPLTEXT
1 0 82 "[3],args[1]+size,args[3],col),Line(args[1]-size,args[4],args[1
]+size,args[4],col);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 4 "end:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 41 "Pex:=Expoint(1.5,2.2,1.9,2.4
,.1,magenta):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "display(\{
Pex\},axes=boxed,view=[-1..2,0..4]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 100 "If desired, one can fill in the circle with more than one, as \+
done here, but it will slow plotting.." }}{PARA 0 "" 0 "" {TEXT -1
335 "We can now plot a circle that when reduced in size will appear f
illed, and an 'error' bar that indicates the uncertainty of the measu
rement. It is straightforward to write routines for other symbols (sq
uares and diamonds at least) We generate a colour table and have the
luxury of distinguishing subsequent data points by colour:" }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 61 "coltable:=[violet,magenta,blue,cyan,green
,yellow,red,maroon];" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "imax
:=7;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 57 "Now we set up a loop to p
repare the data point structure:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "
P3:=[]:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "for i from 1 to imax do:
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "dline:=data[i]:" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 81 "P3:=[Expoint(dline[1],dline[2],dline[3],dline[
4],0.01,coltable[i+1]),op(P3)]: od:" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 23 "display(P3,axes=boxed);" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 99 "With the pro
cedures defined in the previous section we can plot the data set read \+
in from the file " }{TEXT 19 5 "draw1" }{TEXT -1 1 "." }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "display(P3,axes=boxed);" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 66 "The next step is to obtain a linear least
squares fit to the data." }}}{SECT 0 {PARA 3 "" 0 "" {TEXT 270 61 "Us
e Maple's stats package to do the linear least squares fit:" }}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 146 "Suppose we wish to fit the data with a l
inear least squares method: Maple has a package for that, but it requ
ires the data in two lists for the " }{TEXT 271 1 "x" }{TEXT -1 5 " a
nd " }{TEXT 272 1 "y" }{TEXT -1 8 " values." }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 12 "with(stats):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 24 "for i from 1 to imax do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 45 "x
val[i]:=data[i][1]: yval[i]:=data[i][2]: od:" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 106 "We generate the solution after some trial and error: M
aple requires a list of lists in the fit procedure!" }}{PARA 0 "> "
0 "" {MPLTEXT 1 0 53 "lxval:=convert(xval,list); lyval:=convert(yval,l
ist);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "sol:=fit[leastsqua
re[[x,y]]]([lxval,lyval]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
12 "assign(sol);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 2 "y;" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 56 "Pfit:=plot(y,x=xval[1]-0.1..
xval[imax]+0.1,color=brown):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 34 "display(\{Pfit,op(P3)\},axes=boxed);" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 98 "A good question: given the uncertainty in the data, what \+
is the uncertainty in the slope and the " }{TEXT 273 1 "y" }{TEXT -1
63 " intercept? How well are the data fitted by the straight line?" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "describe[linearcorrelation](lxval,
lyval): evalf(%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 157 "The linear \+
correlation is strong (the number is close to 1; 0 would mean no line
ar relationship; -1 that an anticorrelation to a linear relationship e
xists." }}{PARA 0 "" 0 "" {TEXT -1 459 "This illustrates nicely that t
he question of fitting a straight line to measured data points is non-
trivial. Note that the uncertainties have not been used in the calcul
ation. One can weigh the data points, if one trusts some of them more \+
than others. For example: the points on the right apparently follow a \+
straight line rather well. The fit pays attention to this fact, but i
t does not 'know' that the two rightmost points have the biggest erro
r bars. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 77 "What assumptions enter the standar
d linear least squares fit? (J.R. Taylor, " }{TEXT 274 30 "Introducti
on to Error Analysis" }{TEXT -1 16 "; chapter 8.2) " }}{PARA 0 "" 0 "
" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 7 "1) the " }{TEXT 275 1 "x
" }{TEXT -1 69 " values are assumed to be very precise (negligible un
certainties); " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 ""
{TEXT -1 7 "2) the " }{TEXT 276 1 "y" }{TEXT -1 157 " values are assum
ed to have equal uncertainties (if not then weighting should be appli
ed) described by a Gaussian with a certain width common to all y_i; \+
" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 421 "3) i
f the data were exact, they would satisfy y_i = A + B x_i , where A,B
are physical constants. If the uncertainties for each measurement y_
i are governed by a Gaussian width sigma , the probability of obatain
ing the entire data set for the physical constants A,B is given as th
e product of probabilities for each point. The product of exponential
s can be expressed as a single exponential with argument -chi^2/2 ." }
{MPLTEXT 1 0 0 "" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT 277 59 "Derivation \+
expressions for the LLSQ fit with uncertainties:" }}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 51 "y:='y': x:='x': i:='i': chi:='chi': A:='A': B:
='B':" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 145 "The expression for the \+
total error is defined according to the data model (linear relationshi
p, equal Gaussian errors in y_i) described above as:" }}}{EXCHG {PARA
0 "> " 0 "" {MPLTEXT 1 0 48 "E1:=chi^2=Sum((y[i]-A-B*x[i])^2/sigma^2,i
=1..N);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 171 "The best values for A
and B are those for which the probability P_\{A,B\}(y_1,...,y_N) is \+
a maximum, or chi-squared a minimum. They are found from the necessar
y condition:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "E2a:=diff(rhs(E1),A
)=0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "E2b:=diff(rhs(E1),B
)=0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "sol:=solve(\{E2a,E2
b\},\{A,B\}):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "assign(sol
);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "simplify(A);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "simplify(B);" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 154 "Everything goes as per textbook, except \+
that the sum of 1 from 1 to N is not simplified to N. This is due to \+
the fact that the unevaluated Sum was used." }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 9 "value(A);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 9
"value(B);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 26 "A1:=subs(N=im
ax,value(A)):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "i:='i': fo
r j from 1 to imax do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "A1:=subs(
x[j]=xval[j],y[j]=yval[j],A1); od:" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 13 "simplify(A1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT
1 0 26 "B1:=subs(N=imax,value(B)):" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 32 "i:='i': for j from 1 to imax do:" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 43 "B1:=subs(x[j]=xval[j],y[j]=yval[j],B1); od:" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "simplify(B1);" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 330 "Thus, we know how the simple linear leas
t squares result was obtained. Now we are interested in understanding
the estimate of uncertainties. The best estimate for sigma is obtain
ed in analogy with the best estimate for A and B: We find the maximum
probability to obtain the measured values \{y_i\} from the necessary
condition. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT
-1 139 "First we indicate schematically how the probability depends o
n sigma: psi is the same as chi except for the sigma^2 that has been r
emoved:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 38 "Prob:=exp(-psi^2/(2*sigm
a^2))/sigma^N;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "E3:=diff(
Prob,sigma)=0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "solsig:=s
olve(E3,sigma);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "E4:=sigm
a^2=solsig[1]^2;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 246 "Thus, the su
m of squares (psi^2) per degree of freedom (N) provides an estimate f
or the uncertainty of the data. In fact, we need to correct for the f
act that A and B are not exactly known, but rather deduced from the d
ata to divide by (N-2). " }}{PARA 0 "" 0 "" {TEXT -1 24 "For our examp
le we have:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "i:='i':" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "sigsq:=evalf(Sum((yval[i]-A1-B1*xva
l[i])^2,i=1..imax)/(imax-2));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 327
"This number represents a mean square deviation between the data point
s and the line fit. That this is a realistic estimate follows from t
aking a typical deviation from the graph and squaring it (e.g., 0.02 \+
for the green point #4). We are, however, more interested in the unce
rtainty in the extracted coefficients A and B. " }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 31 "evalf(yval[4]-A1-B1*xval[4])^2;" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 226 "Now that we have an error estimate for the data \+
\{y_i\} we can use standard error propagation to estimate the errors \+
in A and B from their respective expressions. Straightforward calculu
s results in the following expressions:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 7 "i:='i':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "
sigBsq:=evalf(sigsq*imax/(imax*Sum(xval[i]^2,i=1..imax)-Sum(xval[i],i=
1..imax)^2));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "i:='i':" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 103 "sigAsq:=evalf(sigsq*Sum(xv
al[i]^2,i=1..imax)/(imax*Sum(xval[i]^2,i=1..imax)-Sum(xval[i],i=1..ima
x)^2));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 120 "These expressions are
not easily derived from A and B. The following example shows Maple's
difficulty with the problem:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "di
ff(sum(z[i]*eta[i],i=1..6),z[2]);" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 35 "diff(sum(z[i]*eta[i],i=1..N),z[2]);" }}}}{SECT 0
{PARA 3 "" 0 "" {TEXT 278 51 "Hand-written LLSQ fit procedure with unc
ertainties:" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 75 "We define procedure
s which implement the derived formulae for the LLSQ fit:" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "addme:=arg->evalf(Sum(arg[i],i=1..n
ops(arg)));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "addme([1,2,3
,4]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "addpr:=(arg1,arg2)
->evalf(Sum(arg1[i]*arg2[i],i=1..nops(arg1)));" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 27 "addpr([1,2,3,4],[1,2,3,4]);" }}}{EXCHG {PARA
0 "> " 0 "" {MPLTEXT 1 0 23 "evalf(sum(i^2,i=1..4));" }}}{EXCHG {PARA
0 "" 0 "" {TEXT -1 68 "We are ready to program the standard linear lea
st squares algorithm:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "i:='i':" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "LLSQ:=proc(xv,yv) local _i,
N,a,b,c,d,Delta,A,B,sq,sigsq,dA,dB;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
12 "N:=nops(xv);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "if nops(yv)<N t
hen RETURN(`Mismatch in sizes of xv and yv`); fi;" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 16 "c:=addpr(xv,xv);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
13 "d:=addme(xv);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "a:=addpr(xv,yv
);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "b:=addme(yv);" }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 22 "Delta:=evalf(N*c-d^2);" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 19 "A:=(c*b-d*a)/Delta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 19 "B:=(N*a-d*b)/Delta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "sq:=[]
; for _i from 1 to N do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "sq:=[op
(sq),evalf(yv[i]-A-B*xv[i])]: od:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
26 "sigsq:=addpr(sq,sq)/(N-2);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "d
A:=sigsq*c/Delta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "dB:=evalf(N*si
gsq)/Delta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "print(`[A,B,sigma,dA
,dB]`);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "[A,B,evalf(sqrt(sigsq),4
),evalf(sqrt(dA),4),evalf(sqrt(dB),4)];" }}{PARA 0 "> " 0 "" {MPLTEXT
1 0 4 "end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "LLSQ([1,2,3,
4],[2.1,4.,5.9,8.2]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "so
l1:=LLSQ(lxval,lyval);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "y
fit:=sol1[1]+sol1[2]*x;" }}}}{SECT 0 {PARA 3 "" 0 "" {TEXT 279 28 "Wei
ghted LLSQ fit procedure:" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 431 "We k
now that our procedure works and we have estimates for the uncertaint
ies in A and B now. We can modify the results to take into account in
dividual uncertainties for the datapoints y_i by means of a weighted \+
least squares fit (J.R. Taylor problem 8.4). The individual uncertain
ties have to be known for this method to be meaningful. One uses for \+
the inverse of the individual sigma_i-squared as a weight for data po
int i." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "i:='i': for i from 1 to \+
imax do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "wval[i]:=1./(0.5*(data[
i][3]-data[i][4]))^2; od:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
26 "lwval:=convert(wval,list);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT
1 0 63 "addtr:=(a1,a2,a3)->evalf(Sum(a1[i]*a2[i]*a3[i],i=1..nops(a1)))
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "WLLSQ:=proc(xv,yv,wv) \+
local Delta,N,a,b,c,d,e,A,B;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "N:=
nops(xv);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 64 "if nops(yv)<N then RET
URN(`Mismatch in sizes of xv and yv`); fi;" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 19 "a:=addtr(wv,xv,xv);" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 16 "b:=addpr(wv,yv);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "c:=addpr(
wv,xv);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "d:=addtr(wv,xv,yv);" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "e:=addme(wv);" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 22 "Delta:=evalf(e*a-c^2);" }}{PARA 0 "> " 0 "" {MPLTEXT
1 0 19 "A:=(a*b-c*d)/Delta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "B:=(
e*d-c*b)/Delta;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "print(`[A,B]`);
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "[A,B];" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 4 "end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "i:='
i':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "WLLSQ(lxval,lyval,lw
val);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 19 "ywlsq:=%[1]+%[2]*x
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "Pwfit:=plot(ywlsq,x=xv
al[1]-0.1..xval[imax]+0.1,color=red):" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 40 "display(\{Pfit,Pwfit,op(P3)\},axes=boxed);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 100 "We see that the data points with \+
larger uncertainties on the right have less control over the slope." }
}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "0 0 0" 0 }
{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }