% EE365: Stochastic Control % simulation of epidemic model for ee365 Markov lecture % Written by Arezou Keshavarz clear all; close all; % this simulation uses the rand() function, and so will give % a different result each time it is run. % to make it repeatable, you can set the rand seed, using rand('state',0) % individuals are located on an n x n grid n = 30; % grid size T = 100; % time steps % A specifies the state of the indviduals. A(i,j,t) is the state of the % individual at the (i,j) location at time t. % 1 is S (susceptible), 2 is I (infected) and 3 is R (removed) A = ones(n,n,T+1); % all individuals susceptible % make the individuals at the (2,2) and (n-1,n-1) location infected. A(2,2,1) = 2; A(n-1,n-1,1) = 2; % for plotting infection spread % S is green, I is red, and R is black. col = ['g','r','k']; cmp = [0 1 0; 1 0 0; 0 0 0]; colormap(cmp) % plot infection spread in population shading faceted pcolor(1:n+1,1:n+1,[A(:,:,1) [1 ;2 ;3; ones(n-3,1)]; ones(1,n+1)]) axis off drawnow; % transition probability matrices % with infected neighbors: % probability of infection is 0.4 if there are infected neighbors % ordering of states is 1,2,3 (S,I,R) P1 = [0.6 0.4 0; 0.2 0.6 0.2; 0 0 1]; % no infected nieghbors % in this case, no chance of infection P2 = [1 0 0 ; 0.2 0.6 0.2; 0 0 1]; %% simulate the infection model for T time steps for t = 1 : T disp(t) for i = 1 : n for j = 1 : n if (A(i,j,t) == 1) % if susceptible, check to see if any neighbors are infected neighbors_I = 0; if (i > 1) neighbors_I = neighbors_I + (A(i-1,j,t) == 2); end if i < n neighbors_I = neighbors_I + (A(i+1,j,t) == 2); end if (j > 1) neighbors_I = neighbors_I + (A(i,j-1,t) == 2); end if (j < n) neighbors_I = neighbors_I + (A(i,j+1,t) == 2); end if neighbors_I > 0 P = P1; else P = P2; end else % infected or removed P = P1; end % generate a sample for the next state, given transiiton % probabilities P(A(i,j,t),:) A(i,j,t+1) = find((cumsum(P(A(i,j,t),:)) - rand) > 0,1); end end % plotting code pcolor(1:n+1,1:n+1,[A(:,:,t+1) [1 ;2 ;3; ones(n-3,1)]; ones(1,n+1)]) drawnow; end %% % generate area plots, to show the percentage of population in the S/I/R state over % time for t = 1 : T+1 total_S(t) = sum(sum(A(:,:,t) == 1))/n^2; total_I(t) = sum(sum(A(:,:,t) == 2))/n^2; total_R(t) = sum(sum(A(:,:,t) == 3))/n^2; end figure(2) plot2 = area(1:T+1,[total_S; total_I; total_R]'); set(plot2(1),'FaceColor',col(1)); set(plot2(2),'FaceColor',col(2)); set(plot2(3),'FaceColor',col(3)); set(gca, 'YLim', [0 1]);