{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 Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 256 "" 1 16 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 
257 "" 0 1 0 0 0 0 1 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 0 0 1 0 0 0 0 0 0 1 }{CSTYLE "" -1 263 "" 0 1 0 0 
0 0 0 0 1 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 0 0 1 0 0 0 
0 0 0 1 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 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 }}
{SECT 0 {EXCHG {PARA 0 "" 0 "" {TEXT 256 21 "Fitting of noisy data" }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 227 "We used \+
a linear least squares fit to analyze the linear relationship between \+
given data. Linear data fitting is applicable for many problems, and i
n fact can also be used for exponential laws (simply take the logarith
m of the " }{TEXT 259 2 "y-" }{TEXT -1 69 "values and then fit), and e
ven power laws (take the logarithm of the " }{TEXT 257 1 "x" }{TEXT 
-1 5 " and " }{TEXT 258 1 "y" }{TEXT -1 190 " values and fit). In this
 worksheet we continue along the same lines, but extend the fit to all
ow for a quadratic function. The fit is still linear in the parameters
 to be determined, i.e., " }{XPPEDIT 18 0 "y = a+bx+cx^2;" "6#/%\"yG,(
%\"aG\"\"\"%#bxGF'*$%#cxG\"\"#F'" }{TEXT -1 404 ". Note that the probl
em where nonlinear parameters are to be determined, and where the prob
lem cannot be reduced to a linear fit of a related problem (as mention
ed above), is considerably more difficult. No built-in procedures are \+
available in Maple for nonlinear fitting at the present time. There ex
ist implementations of some algorithms for this case described in the \+
literature (cf. W. Press et al., " }{TEXT 270 17 "Numerical Recipes" }
{TEXT -1 30 ", Cambridge University Press)." }}{PARA 0 "" 0 "" {TEXT 
-1 408 "We follow another approach in the following respect: we will u
se a random number generator to produce data with uncertainties of a k
nown variance. Then we will test the fitting techniques and the error \+
estimates. Thus, we will use Maple to perform a computer experiment. W
e will rely solely on Maple's built-in fit procedure, even though it d
oes not provide us with uncertainties for the fitted coefficients." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 136 "We start
 with the drawing tools that we developed for plotting symbols with er
ror bars. They work reasonably well if the scaling on the " }{TEXT 
261 1 "x" }{TEXT -1 5 " and " }{TEXT 260 1 "y" }{TEXT -1 143 " axes is
 not very different. You may need to scale the data if they don't perf
orm well for you (or replace the Circle procedure by an Ellipse)." }}}
{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 112 "plot([args[3]*co
s(t)+args[1],args[3]*sin(t)+args[2],t=0..2*Pi],scaling=CONSTRAINED,col
our=col,thickness=2); end:" }}}{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 80 "plot([[args[1],args[2]],[args[3],args[4]]],scaling=CO
NSTRAINED,colour=col); end:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 
48 "Expoint:=proc() local size,col; size:=2*args[5];" }}{PARA 0 "> " 
0 "" {MPLTEXT 1 0 48 "if nargs<6 then col:=black else col:=args[6] fi;
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 222 "Circle(args[1],args[2],0.5*arg
s[5],col),Circle(args[1],args[2],args[5],col),Line(args[1],args[3],arg
s[1],args[4],col),Line(args[1]-size,args[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 
12 "with(plots):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 236 "We need a ra
ndom number generator that produces random numbers on the interval [-1
,1], which are then easily transformed to some other interval. Maple h
as by default an integer random number generator. The stats package is
 of some help." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(stat
s);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "with(stats[statplots
]):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 366 "We learn that there is a \+
random number generator for normally or uniformly distributed data, an
d a histogram function that works either with tallied (adjusted) inter
vals so that each bin has equal area, or for equidistant bins. The bin
 size is calculated by the program, and the area under the histogram e
quals unity to allow for a direct probability interpretation." }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "histogram([random[normald](2
50)]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "data:=random[norm
ald](250):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "histogram[1](
[data]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "xrand:=stats[ra
ndom, uniform](20);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "hist
ogram[1]([xrand]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 57 "These are (
supposedly) uniformly distributed over [0,1]. " }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 262 16 "Linear fit case:" }}
{PARA 0 "" 0 "" {TEXT -1 58 "We generate 5 data points with random err
ors according to " }{TEXT 265 2 "y " }{TEXT -1 10 "= 1 + 0.5 " }{TEXT 
264 1 "x" }{TEXT -1 13 ". We use the " }{TEXT 266 1 "x" }{TEXT -1 15 "
 interval 0..2:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "Y:=x->1+
0.5*x;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "errsize:=0.1;" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "imax:=5: for i from 1 to im
ax do: xval[i]:=2/imax*(i-0.5):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "
yval[i]:=Y(xval[i])+errsize*(2.*xrand[i]-1.); od:" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 32 "P1:=plot(Y(x),x=0..2,color=red):" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "P2:=\{\}: i:='i': symbolsize:=0.02:
 for i from 1 to imax do:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 102 "P2:=
\{op(P2),Expoint(xval[i],yval[i],yval[i]+errsize,yval[i]-errsize,symbo
lsize,blue)\}: od: P2:=op(P2):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 
1 0 17 "display(\{P1,P2\});" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 184 "N
ote that our error bounds are not standard errors, but brackets on the
 error as we have a uniform distribution of the random numbers used to
 produce the data points from exact values." }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 214 "We are ready to perform a leas
t squares fit. We can get help by requesting ?fit. We need to do some \+
experimenting to realize that fit cannot take tables, but needs lists \+
of data values. Thus we convert in the call:" }}{PARA 0 "" 0 "" {TEXT 
-1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "y:='y': sol:=fit[
leastsquare[[x,y], y=a*x+b]]( [convert(xval,list),convert(yval,list)])
;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "assign(sol);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "P3:=plot(y,x=0..2,color=gree
n):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "display(\{P1,P2,P3\}
);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 15 "For comparison:" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "sort(y); sort(Y(x));" }}}{EXCHG 
{PARA 0 "" 0 "" {TEXT -1 32 "The difference in the slope and " }{TEXT 
267 1 "y" }{TEXT -1 131 " intercept between the 'true' and fitted resu
lts is small, in fact much smaller than the absolute error bound for e
ach data point. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT 263 19 "Quadratic fit case:" }}{PARA 0 "" 0 "" {TEXT -1 
0 "" }}{PARA 0 "" 0 "" {TEXT -1 58 "We generate 6 data points with ran
dom errors according to " }{XPPEDIT 18 0 "y = 1-.5*x+.5*x^2;" "6#/%\"y
G,(\"\"\"F&*&$\"\"&!\"\"F&%\"xGF&F**&$F)F*F&*$F+\"\"#F&F&" }{TEXT -1 
2 ". " }}{PARA 0 "" 0 "" {TEXT -1 11 "We use the " }{TEXT 268 1 "x" }
{TEXT -1 139 " interval 0..2 again, and keep the error size. The numbe
r of data points was increased, since there is an extra parameter to b
e determined." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 22 "Y:=x->1-0.
5*x+0.5*x^2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "imax:=6: fo
r i from 1 to imax do: xval[i]:=2/imax*(i-0.5):" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 49 "yval[i]:=Y(xval[i])+errsize*(2.*xrand[i]-1.); od:" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "P1:=plot(Y(x),x=0..2,color=
red):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "P2:=\{\}: i:='i': \+
symbolsize:=0.02: for i from 1 to imax do:" }}{PARA 0 "> " 0 "" 
{MPLTEXT 1 0 102 "P2:=\{op(P2),Expoint(xval[i],yval[i],yval[i]+errsize
,yval[i]-errsize,symbolsize,blue)\}: od: P2:=op(P2):" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 17 "display(\{P1,P2\});" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 184 "Note that our error bounds are not standard error
s, but brackets on the error as we have a uniform distribution of the \+
random numbers used to produce the data points from exact values." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 100 "We are r
eady to perform a least squares fit. The model is now specified to inc
lude a quadratic term:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 94 "y:='y': sol:=fit[leastsquare[[x,y],
 y=a*x^2+b*x+c]]( [convert(xval,list),convert(yval,list)]);" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "assign(sol);" }}}{EXCHG 
{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "P3:=plot(y,x=0..2,color=green):" }}
}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "display(\{P1,P2,P3\});" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 15 "For comparison:" }}}{EXCHG {PARA 
0 "> " 0 "" {MPLTEXT 1 0 20 "sort(y); sort(Y(x));" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 210 "We learn that a least squares fit for a quadratic
 function can be a trickier problem: the added flexibility puts a grea
ter demand on the problem. The coefficients deviate more from the 'tru
e' model (red curve)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 
"" {TEXT 269 10 "Exercises:" }{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT 
-1 142 "1) Improve the situation by adding more data points in the int
erval. Also try fewer data points. Repeat the fits with increased erro
r margins." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 
-1 178 "2) Observe the deviation between the 'true' and the model resu
lt based on randomized data fitting outside the interval in which the \+
fitting was performed (extrapolation problem)." }}{PARA 0 "" 0 "" 
{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 127 "3) Repeat the data fitti
ng problems with normally distributed errors (the error size plays the
 role of the standard deviation)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 50 "4) Test the data fit problem for a cubic \+
function." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{MARK "25
 0 0" 20 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 
}