Page 60
mcent:=[3.6, 3.7, 3.9, 4.3, 4.2, 4.0, 3.5, 2.9,
   2.2,2.0, 2.2, 2.1, 1.9, 1.8,1.4, 0.8, 0.3]:
fcent:=[3.6, 3.9, 4.1, 4.7, 5.0, 4.3, 4.0, 3.6,
   2.7, 2.8, 3.0, 3.1, 2.8, 2.3, 2.0, 1.7, 1.6]:
tot:=[seq(mcent[i]+fcent[i],i=1..17)]:
ranges:=[0..5, 5..10, 10..15, 15..20,
   20..25, 25..30, 30..35, 35..40, 40..45,
   45..50, 50..55, 55..60, 60..65, 65..70,
   70..75, 75..80, 80..85]:
with(stats): with(plots):
mpop:=[seq(Weight(ranges[i], 5*mcent[i]),i=1..17)]:
fpop:=[seq(Weight(ranges[i], 5*fcent[i]),i=1..17)]:
pop:=[seq(Weight(ranges[i], 5*tot[i]),i=1..17)]:
statplots[histogram](pop);

              Page 61
age:=[2.5, 7.5, 12.5, 17.5, 22.5, 27.5, 32.5,37.5
   42.5, 47.5, 52.5, 57.5, 62.5, 67.5, 72.5, 77.5, 82.5];
cummale:=[seq(sum('mcent[i]','i'=1..n),n=1..17)]:
cumfale:=[seq(sum('fcent[i]','i'=1..n),n=1..17)]:
cumtot:=[seq(sum('tot[i]','i'=1..n),n=1..17)]:
ptsm:=[seq([age[i],cummale[i]],i=1..17)];
ptsf:=[seq([age[i],cumfale[i]],i=1..17)]:
ptsT:=[seq([age[i],cumtot[i]],i=1..17)]:
plot({ ptsm,ptsf,ptsT},color=BLACK);

              Page 63
Sum('age[j]'*'tot[j]',j=1..17)= sum(age[j]*tot[j],j=1..17);
Sum('mcent[n]*age[n]','n'=1..17)/Sum('mcent[n]','n'=1..17)
  = sum('mcent[n]*age[n]','n'=1..17)/sum('mcent[n]','n'=1..17);
with(describe): mean(pop); median(pop);

              Page 64
bus1:=[60,60,60,60]; bus2:=[30,0,120,90];
range(bus1), range(bus2);
median(bus1), median(bus2);
mean(bus1), mean(bus2);
variance(bus1), variance(bus2);
standarddeviation(bus1), standarddeviation(bus2);

range(mpop), range(fpop), range(pop);
median(mpop), median(fpop), median(pop);
mean(mpop), mean(fpop), mean(pop);
variance(mpop), variance(fpop), variance(pop);
standarddeviation(mpop), standarddeviation(fpop),
standarddeviation(pop);

             Page 65
y:=(sigma,mu,x)-> exp(-(x-mu)^2/
   (2*sigma^2))/(sqrt(2*Pi)*sigma);
plot({ y(1,0,x),y(2,0,x),y(3,0,x)},x=-10..10);
plot({ y(1,-4,x),y(1,0,x),y(1,4,x)} ,x=-10..10);

              Exercises/Experiments

1.
with(stats): with (plots): with(describe):
Mort:=[1122.4, 55.1, 27.5, 33.4, 118.4, 139.6, 158.0, 196.4,
 231.0, 287.8, 487.2, 711.2, 1116.9, 1685.1, 2435.5, 3632.4,
 5300.0, 8142.0, 15278.0]:
MortRate:=[seq(Mort[i]/100000,i=1..19)];
ranges=[seq(5*i..5*(i+1),i=1..17)];
mortdata=[Weight(0..1,MortRate[1]), Weight(1..5,4*MortRate[2]),
  seq(Weight(ranges[i],5*MortRate[2+i]), i=1..17)]:
statplots[histogram](mortdata);

1a.
xcord:=[seq(3+5*(i-1),i=1..18)];
mortrate:=[seq(MortRate[i+1], i=1..18)];
plot([seq([xcord[i],mortrate[i]],i=1..18)],
   style=POINT, symbol=CROSS);
fit[leastsquare[[x,y],y=a+b*x+c*x^2+d*x^3]]([xcord,mortrate]);
approx:=unapply(rhs(%),x);
approx(30)*100000;
plot(approx(x),x=0..90);

1b.
Lnmortrate:=map(ln,mortate);
fit[leastsquare[[x,y],y=m*x+b]]([xcord,Lnmortrate]);
k:=op(1,op(1,rhs(%))); A:=op(2,rhs(%%));
expfit:=t-> exp(A)*exp(k*t);
expfit(30)*100000;
J:=plot(expfit(t),t=0..85):
K:=plot([seq([xcord[i],MortRate[i+1]],i=1..18)],
  style=POINT,symbol=CROSS):
display({ J,K});

1c.
readlib(spline):
linefit:=spline(xcord,mortrate,x,linear):
y:=unapply(linefit.x): y(30)*100000;
J:=plot(y(t), t=0..85):
display({ J,K});

summary:=[Weight(Mort[1],1),Weight(Mort[2],4),
  seq(Weight(Mort[i],5),i=3..19)];
range(summary); median(summary); mean(summary);
  standarddeviation(summary);

2.
with(stats): with(plots): with(describe):
htinches:=[seq(60+i,i=1..15)];
numMales:= [0,0,0,0,0,2,1,2,7,10,14,7,5,2,1];
ranges:=[seq(htinches[i]..htinches[i]+1,
   i=1..15)];
maledata:=[seq(Weight(ranges[i],numMales[i]),
   i=1..15)];
statplots[histogram](maledata);
range(maledata); median(maledata);
 mean(maledata); standarddeviation(maledata);
'The average height is',floor(%%/12),
'feet and',floor(frac(%%%/12)*12),
'inches';
'The standard deviation is',
   floor(frac(%%%/12)*12), 'inches';

mu:=mean(maledata);
sigma:=standarddeviation(maledata);
ND:=x-> exp(-(x-mu)^2/(2*sigma^2))/(sigma*sqrt(2*Pi));
J:=plot(mu*ND(x),x=60..76):
K:=statplots[histogram](maledata):
plots[display]({ J,K});

4.
sigma:='sigma ': mu='mu ':
f:=x-> exp(-(x-mu)^2/(2*sigma^2))/(sigma*sqrt(2*Pi));
assume(sigma > 0);
int(x*f(x),x=-infinity..infinity);
int((x-mu)^2*f(x),x=-infinity..infinity);