%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % % A solution to CSE Your Homework Assignment Project 5 % % Models of Infection: Person to Person % % Dianne P. O'Leary 09/03 % % % % solution.m % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % This project investigates a model for the spread of % % infection in a hospital ward. % % % % The parameters for the model: % % % % n,m: dimensions of the grid of hospital beds % % k: number of days of contagion % % tau: transmission rate for the infection % % delta: mobility rate of the population % % nu: vaccination rate % % % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % The following global variables are passed to % % function zeroinfect, for Problem 5. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% global target ntimes n k tau delta movieon = 1; movieoff = 0; % Plotting symbols: pl(1).s = 'b'; pl(2).s = 'r-.'; pl(3).s = 'k:'; pl(4).s = 'g--'; pl(5).s = 'm'; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Problem 1 % % Run a Monte Carlo simulation with no mobility or % % vaccination. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(1) n = 10; mn = n*n; k = 4; tau = .2; delta = 0; nu = 0; [nrecovered1,infected,recovered,susceptible,vaccinated] = ... modelmc(n,k,tau,delta,nu,movieon); t = length(infected); plot(1:t,infected/mn,pl(2).s,1:t,susceptible/mn,pl(3).s, ... 1:t,recovered/mn,pl(1).s,'LineWidth',1.5) legend('Infected','Susceptible','Recovered') xlabel('day') ylabel('Proportion of individuals') title(sprintf('Disease Status with tau = %f',tau)) print -depsc figure1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Problem 2 % % Add mobility of the patients to the model. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(2) delta = .01; [nrecovered2,infected,recovered,susceptible,vaccinated] = ... modelmc(n,k,tau,delta,nu,movieoff); t = length(infected); plot(1:t,infected/mn,pl(2).s,1:t,susceptible/mn,pl(3).s, ... 1:t,recovered/mn,pl(1).s,'LineWidth',1.5) legend('Infected','Susceptible','Recovered') xlabel('day') ylabel('Proportion of individuals') title(sprintf('Disease Status with tau = %f, delta = %f' ... ,tau, delta)) print -depsc figure2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Problem 3 % % Add vaccination to the model. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(3) nu = .1; [nrecovered3,infected,recovered,susceptible,vaccinated] = ... modelmc(n,k,tau,delta,nu,movieoff); t = length(infected); plot(1:t,infected/mn,pl(2).s,1:t,susceptible/mn,pl(3).s, ... 1:t,recovered/mn,pl(1).s,1:t,vaccinated/mn,pl(4).s,'LineWidth',1.5) legend('Infected','Susceptible','Recovered','Vaccinated') xlabel('day') ylabel('Proportion of individuals') title(sprintf('Disease Status with tau = %f, delta = %f, nu= %f' ... ,tau, delta,nu)) print -depsc figure3 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Problem 4 % % Investigate the variance of the estimate of % % the number of patients eventually infected % % (nrecovered) by taking 1000 trials using various % % nu values. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(4) jj = 1; for nu=0:.1:.3, for i=1:1000, recovered(i) = modelmc(n,k,tau,delta,nu,movieoff)/mn; end subplot(2,2,jj) hist(recovered*100) xlabel('Percent infected') ylabel('Number of trials') title(sprintf('Histogram of infection rate for nu = %f',nu)) disp(sprintf('For nu= %f, Mean = %f, Variance = %f', ... nu,mean(recovered),var(recovered))) jj = jj + 1; end print -depsc figure4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Problem 5 % % Compute the vaccination rate that confines the % % expected infection level to 20 per cent of the % % population. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% figure(5) target = .2; ntimes = 100; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % We can try to use the zerofinder fzero, but the % % results are unreliable. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% nu = fzero(@zeroinfect,.1) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Problem 6 % % Investigate a Markov Chain model for a 2-day % % infection on a 3-person population. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% data = [ 2 1 (1-tau)^2 3 1 tau*tau 5 1 tau*(1-tau) 8 1 tau*(1-tau) 13 2 tau*(1-tau) 14 2 tau*tau 15 2 tau*(1-tau) 17 2 (1-tau)^2 4 3 1 18 4 1 6 5 tau 11 5 1-tau 7 6 1 18 7 1 9 8 tau 12 8 1-tau 10 9 1 18 10 1 19 11 1 16 12 1 12 13 1 4 14 1 11 15 1 16 16 1 17 17 1 18 18 1 19 19 1]; a = sparse(data(:,1),data(:,2),data(:,3)); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compute a stationary vector for the matrix. % % Since the matrix is reducible, there are many % % stationary vectors, so we isolate the one of interest % % by using the power method started with a vector z % % corresponding to starting in the first state. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% z = zeros(19,1); znew = z; znew(1)=1; while norm(z-znew) > eps, z = znew; znew = a*z; end prob(1) = znew(17); prob(2) = znew(16)+znew(19); prob(3) = znew(18); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compare with results of Monte Carlo trials. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% i=0; probmc = zeros(1,3); while (norm(probmc-prob,'inf')>.005) i=i+1; nrecovered(i) = modelmc1d(3,2,tau); probmc(1) = sum(nrecovered==1)/i; probmc(2) = sum(nrecovered==2)/i; probmc(3) = sum(nrecovered==3)/i; end ntrials = i; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Compare with analytical results. % %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% probtr(1) = (1-tau)^4; probtr(2) = 2*(tau*(1-tau)^2 + tau*(1-tau)^3); probtr(3) = tau^2 + 2*tau^2*(1-tau) + (1-tau)^2*tau^2; disp('The 3-bed case') for i=1:3, disp(sprintf('probability of %d infection predicted to be',i)) disp(sprintf( ... ' %f by Markov chain, %f by Monte Carlo, and %f analytically',... prob(i),probmc(i),probtr(i))) end disp(sprintf(... 'It took %d trials of Monte Carlo to get 2 digits of accuracy',... ntrials))