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];
total=mcent+fcent;
x=[5:5:85];  % 5, 10, 15, ..., 85
bar(x,total) % bars centered on the x values
xlabel('Age(years)')
ylabel('Percent in age bracket');
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];
total=mcent+fcent;
x=[5:5:85];  % 5, 10, 15, ..., 85
bar(x,total) % bars centered on the x values
xlabel('Age(years)')
ylabel('Percent in age bracket');

                   Page 61
cumM=cumsum(mcent);
cumF=cumsum(fcent);
cumTot=cumsum(total)
plot(x,cumTot,x,cumM,x,cumF)

                   Page 63
xmid=[2.5:5:82.5];
pop=xmid.*total; % term by term mult.=percentage wt'd ranges
muTotal=sum(pop)/100 %divide by 100 as data is in percent
muM=sum(xmid.*mcent)/sum(mcent)
muF=sum(xmid.*fcent)/sum(fcent)

                   Page 64
bus1=[60 60 60 60]; bus2=[30 0 120 90];
max(bus1), min(bus1)
max(bus2), min(bus2)
median(bus1), median(bus2)
mean(bus1), mean(bus2)
cov(bus1), cov(bus2)
std(bus1), std(bus2)

v=(xmid-mu).^2
   % unweighted vector of deviations squared
var=sum(v.*total)/100 % variance of total pop.
sqrt(var) % std dev of the total pop.

                   Page 65
% make up an m-file, gaussian.m:
   % function y=gaussian(x,m,s);
   %   % m=mean, s=stddev
   %   % note 1/sqrt(2*pi)=.3989422803
   % y=(.3989422803/s)*
   % exp(-0.5*((x-m)./s).^2);
x=[-10:.1:10];
y=gaussian(x,0,1); plot(x,y);hold on;
y=gaussian(x,0,2); plot(x,y);
y=gaussian(x,0,4); plot(x,y);
hold off
y=gaussian(x,0,1); plot(x,y);hold on
y=gaussian(x,-5,1); plot(x,y);
y=gaussian(x,5,1); plot(x,y);

                   Exercises/Experiments

1.
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=Mort/1000;
x=[.5,2.5:5:87.5];
bar(x,MortRate)
x=x(2:19) % 1st point an outlier
MortRate=MortRate(2:19) % ditto

1a.
% cubic fit rate = d*x^3+c*x^2+b*x+a
p=polyfit(x,MortRate,3)
% use builtin polynomial fitter, 3rd order
y=polyval(p,x); % fit evaluated at the x's
plot(x,MortRate,'x'); hold on
plot(x,y)
% or use the general leastsquares model
MT=[x.^3; x.^2; x; ones(size(x))];
cubic=MT'\MortRate'
y=polyval(cubic,x); plot(x,y)

1b.
% exponential fit log(MortRate)=
%   a+b*x or MortRate=exp(a)*exp(bx)
Lnmortrate=log(MortRate);
MT=[ones(size(x)); x];
expon=MT'\Lnmortrate'
hold off
plot(x,MortRate,'x'); hold on
plot(x,exp(expon(1))*exp(expon(2)*x))

1c.
% linear spline fit = straight line
% between points, usual matlab method
hold off
plot(x,MortRate,x,MortRate,'x')

% to interpolate any desired value, use
% interp1, e.g. rate=interp1(x,MortRate,70)
%   interpolated value at x=70
%  mean, median and std dev (of Mortality
%    weighted by age)
%
size(Mort)
wt=[1,4,5*ones(1,17)]
wtSum = Mort*wt' % dot product
mu=wtSum/sum(wt)
median(Mort) % picks out the middle value, no duplicates here
v=(Mort-mu).^2; % vector of squared differences
var=sum(v.*wt)/sum(wt);
std=sqrt(var)

2.
htinches=61:75;
numMales=[0,0,0,0,0,2,1,2,7,10,14,7,5,2,1];
bar(htinches,numMales)
min(htinches)
max(htinches) % range = from min to max
unrolled=[];  % dup. each ht by its # cases
s=size(htinches);
for k=1:s(2)
   j=numMales(k);
   while j>0
     unrolled=[unrolled, htinches(k)];
     j=j-1;
   end
end
median(unrolled)
mu=mean(unrolled+.5) % eg ht 66 counts as 66.5
% alternatively
mu=dot((htinches+.5),numMales)/sum(numMales)
v=(htinches+.5-mu).^2;
var=sum(v.*numMales)/sum(numMales)
std=sqrt(var)

x=60:.1:76;
y=exp(-((x-mu)/std).^2/2)/(std*sqrt(2*pi));
bar(htinches,numMales/sum(numMales))
hold on; plot(x,y)

4.
% to integrate with matlab one can use
% trapz(x,y) on x and y vectors or use
% Simpson's rule, quad(), but this
% requires an m-file.
%  Here we will use the trapzoidal rule.
x=linspace(-3,3);
% simulates -infinity to +infinity here
y=exp(-x.^2/2)/sqrt(2*pi);
trapz(x,y)  % approximately 1