Page 360
% make an m-file, incubationProfile.m:
% function y=incubationProbile(t);
% if t<2
%   y=0;
% elseif t<18
%   y=1/16;
% else
%   y=0;
% end
t=linspace(0,20,100)
for k=1:100
   p1(k)=incubationProfile(t(k)); % uniform
end
plot(t,p1); hold on
p2=exp(-t/6)/6; plot(t,p2) % exponential
p3=t.^9.*exp(-t);
c1=trapz(t,p3); % for normalization
p3=p3/c1; plot(t,p3) % gamma distribution
for k=1:100
   if t(k)<2
      p4(k)=0;
   elseif t(k)<20
     p4(k)=sqrt((t(k)-2)/16)*(1-(t(k)-2)/16)^4;
   else
      p4(k)=0;
   end
end
c2=trapz(t,p4);
p4=p4/c2; plot(t,p4) % the beta distribution

               Page 364
% make an m-file, casesOnset.m:
% function a=casesOnset(t);
%   if t<0
%     a=0;
%   elseif t<2
%     r=linspace(0,t,20);
%     for i=1:20
%       y(i)=incubationProfile(r(i));
%     end
%     a=1000*trapz(r,y);
%   else
%     r=linspace(t-2,t,20);
%    for i=1:20
%        y(i)=incubationProfile(r(i));
%     end
%     a=1000*trapz(r,y);
%   end  % end of casesOnset.m
% recall t=linspace(0,20,100) 
for k=1:100
   a(k)=casesOnset(t(k))   % vector same size as t
end
plot(t,a)
trapz(t,a) % integral 0 to 20

                 Page 366
% for onset of cases with the gamma incubation period,
% make up an m-file, onsetGam.m, containing
% function a=onsetGam(t);
% c1=3.6107e+05;  % from text calculation
% if t<0
%    a=0;
%  elseif t<2
%    r=linspace(0,t,20);
%   for i=1:20
%     y(i)=r(i)^9*exp(-r(i))/c1; %values of t.^9.*exp(-t) at r(i)
%   end
%   a=1000*trapz(r,y);
% else
%   r=linspace(t-2,t,20); % integral of y from 0 to t
%   for i=1:20
%     y(i)=r(i)^9*exp(-r(i))/c1; %values of t.^9.*exp(-t) at r(i)
%   end
%   a=1000*trapz(r,y);
% end
t=linspace(0,20,100);
for k=1:100
   a(k)=onsetGam(t(k)); % creates a vector same size as t
end
plot(t,a)
trapz(t,a) % integral of a over 0 to 20


                 Exercises/Experiments
2a.
% contents of m-file infectDensity.m:
% function h=infectDensity(t)
% % won't vectorize due to the conditionals
% if t<1980
%    h=0;
% elseif t<2020
%    h=((t-1980)/40)*(1-((t-1980)/40)^6);
% else
%    h=0;
% end
t=linspace(1980,2020); % 100 values
for k=1:100
   h(k)=infectDensity(t(k));
end
plot(t,h)

2b.
hmax=max(h)
m=0;
for i=1:100
  if h(i)==hmax
     m=i;
   end
end
maxYear=t(m)

2c.
s=linspace(2,18);
incDen=((s-2)/16).*(1-((s-2)/16).^2);
c5=trapz(s,incDen)
incDen=incDen/c5;
plot(r,incDen)

2d.
% the integral a(t)=int(h(t-s)*p(s)*ds)
% from 0 to infinity has 5 cases:
%     t<1982, a(t)=0,
% 1982<t<1998, a(t)= int(h(t-s)*p(s)ds from 2 to t-1998
% 1998<t<2022, a(t)= int(h(t-s)*p(s)ds from 2 to 18
% 2022<t,      a(t)= int(h(t-s)*p(s)ds from t-2020 to 18
% Here's why.  First case: suppose t=1981; since p(s)=0 unless s>2,
% t-s<1979 so h(t-s)=0.
% Second case: suppose t=1994;
% again p(s)=0 unless s>2, so the lower limit must be at least 2.
% For t=1994, t-1980=14, so s runs from 2 to 14 and t-s runs from
% 1992 down to 1980; after that h(t-s)=0.
% We leave the remaining cases for you.

% But all cases will be done automatically since we have defined
% infectDensity(t)=0 for t<1980 or t>2020
t=linspace(1980,2050);  % get the graph for 1980 to 2050
for k=1:100  % k= time index
   T=t(k);
   for j=1:100 % j= s index
     S=s(j);
     y(j)=infectDensity(T-S)*incDen(j);
   end
   a(k)=trapz(s,y);
end
plot(t,a)