Page 20
ht:=[75,92,108,121,130,142,155];
age:=[1,3,5,7,9,11,13];
sumy:=sum(ht[n],n=1..7); sumx:=sum(age[n],n=1..7);
sumx2:=sum(age[n]^2,n=1..7);
sumxy:=sum(age[n]*ht[n],n=1..7);
m:=evalf((7*sumxy-sumx*sumy)/(7*sumx2-sumx^2));
b:=evalf((sumx2*sumy-sumx*sumxy)/(7*sumx2-sumx^2));



         Page 21
m:='m'; b:='b';# clears m and b
# next create an array of (age,ht) pairs;
pts:=[seq([age[i],ht[i]],i=1..7)]:
with(plots): with(stats):
Data:=plot(pts,style=POINT,symbol=CIRCLE):
fit[leastsquare[[x,y],y=m*x+b]]([age,ht]);
# result in y=m*x+b form, m*x is the
# 1st operand on the right hand side
m:=op(1,op(1,rhs(%))); # strip off x too
b:=op(2,rhs(%%)); # use %% to get 2nd statement back
Fit:=plot(m*x+b,x=0..14):
display({ Data,Fit});



       Page 22/23
restart:
AIDS:=([97, 206, 406, 700, 1289, 1654, 2576,
 3392, 4922, 6343, 8359, 9968, 12990, 14397,
 16604, 17124, 19585, 19707, 21392, 20846,
 23690, 24610, 26228, 22768, 4903]);
CAC:=[seq(sum(AIDS[j]/1000.0, j=1..i),i=1..24)];
Time:=[seq(1981+(i-1)/2,i=1..24)]:

LnCAC:=map(ln,CAC);
Lntime:=map(ln,[seq((i+1)/2/10,i=1..24)]);

with(stats):
fit[leastsquare[[x,y],y=k*x+LnA]]([Lntime,LnCAC]);
k:=op(1,op(1,rhs(%))); LnA:=(op(2,rhs(%%)));
A:=exp(LnA);

Lndata:=plot([seq([Lntime[i],LnCAC[i]],i=1..24)]
  style=POINT,symbol=CIRCLE):
Lnfit:=plot(k*x+ln(A),x=-2.5..0.5):
plots[display]({ Lndata,Lnfit} );

      Page 24
n:=trunc(k);

pts:=[seq([Time[i], CAC[i]], i=1..24)];
Fit:=plot(A*((t-1980)/10)^n,t=1980..1993):
Data:=plot(pts,style=POINT,symbol=CIRCLE):
plots[display]({Fit,Data});




         Exercises/Experiments

1.
ht:=[62,63,64,65,66,67,68,69,70,71,72,73,74];
wt:=[128,131,135,139,142,146,150,154,158,162,167,172,177];
with(stats):
fit[leastsquare[[x,y], y= a*x^3+b*x^2+c*x+d]]([ht,wt]);
y:=unapply(rhs(%),x);
pts:=[seq([ht[i],wt[i]],i=1..13)];
J:=plot(pts,style=POINT,symbol=CROSS): K:=plot(y(x),x=62..74):
with(plots): display({ J,K});
errorLinear:= sum('(4.04*ht[i]-124.14- wt[i])^2','i'=1..13);
errorcubic:= sum('(y(ht[i])- wt[i])^2','i'=1..13);
evalf(%);

2.
restart:
age60:=[0,10,20,30,40,50,60,80,100]:
percent60:=[100,98.5,98,96.5,95,92.5,79,34,4]:
with(stats):
fit[leastsquare[[x,y],y=a*x^4+
  b*x^3+c*x^2+d*x+e]]([age60,percent60]);
yfit60:=unapply(rhs(\qu),x):
pts60:=[seq([age60[i],percent60[i]],i=1..9)]:
J6:=plot(pts60,style=POINT,symbol=CROSS):
K6:=plot(yfit60(x),x=0..100):
with(plots): display({ J6,K6});

3.
restart:
AIDS:=([97, 206, 406, 700, 1289, 1654,2576, 3392, 4922, 6343, 8359,
 9968, 12990, 14397, 16604, 17124, 19585,19707, 21392, 20846, 23690,
 24610, 26228, 22768, 4903]);
CAC:=[seq(sum(AIDS[j]/1000.0,j=1..i),i=1..24)];
Time:=[seq(1981+(i-1)/2,i=1..24)]:
pts:=[seq([Time[i],CAC[i]],i=1..24)]:
LnCAC:=map(ln,CAC);
Times:=[seq((i+1)/2/10,i=1..24)];
with(stats):
fit[leastsquare[[x,y],y=m*x+b]]([Times,LnCAC]);
k:=op(1,op(1,rhs(%))); A:=op(2,rhs(%%));
y:=t-> exp(A)*exp(k*t);
J:=plot(y((t-1980)/10),t=1980..1992):
K:=plot(pts,style=POINT,symbol=CIRCLE):
plots[display]({J,K});}

4.
CW:=[24.2,22.9,27.,21.5,23.5,22.4,
 23.8, 25.5, 24.5,25.5,22.,24.5];
GSW:=[38.5,26.,34.,25.5,37.,30.,34.,
 43.5,30.5, 36.,29.,32];
with(stats):
fit[leastsquare[[x,y],y=m*x+b]]([CW,GSW]);
pts:=[seq([CW[i],GSW[i]],i=1..12)];
J:=plot(pts,style=POINT,symbol=CROSS):
K:=plot(2.107*x-17.447,x=21..28):
CM:=[28.5,24.5,26.5,28.25,28.2,29.5,
  24.5,26.9,28.2,25.6,28.1,27.8,29.5,29.5,29];
GSM:=[45.8,47.5,50.8,51.5,55.0,51.,47.5,45.,
  56.0,49.5,57.5,51.,59.5, 58.,68.25];
fit[leastsquare[[x,y],y=m*x+b]]([CM,GSM]);
pts:=[seq([CM[i],GSM[i]],i=1..15)];
L:=plot(pts,style=POINT,symbol=CIRCLE):
M:=plot(2.153*x-6.567,x=24..30):
with(plots): display({ J,K,L,M});