{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) " 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) " 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 }