solve using MATLAB. show code EXERSICES PROBLEMS 1.-Write a program on a worksta

Question

solve using MATLAB. show code

EXERSICES PROBLEMS 1.-Write a program on a workstation that emulates the mathematical epidemic model of Reed and Frost. For purposes of simplicity, fix the number of individuals in the population at 1,000 and fix the maximum number of time periods at 25. However, allow the probability of effective contact and the number of immunes to vary s that the effects of social service programs such as immunization or prevention education can be observed. Thus, the only input to your program will be the initial number of immunes and the effective contact probability For output, your program should display in each time period the number of 1. susceptibles 2. active cases 3. immunes Run the program many times while varying the number of immunes and the probability of effective contact. Sensitize yourself to epidemic patterns and variations.

Explanation / Answer

% clear; close all

profile on

%Transmission network generation

k=5;

N=500;

g=1;

tic

% T=make_homogeneous_symmetric_network(N,k);%note that this does not always converge to a

solution

% T=make_random_network(N,k);%

% T=make_lattice_network(M,N);

% T=make_scale_free(N,k); %k is min number contacts per individual

% T=make_small_world;

% Predicted_R0=tau*k/g

No_sims=2000;

outMat=[];

tau=[.31];

no_trials=length(tau);

Frac_take=zeros(2,length(tau));

for i=1:no_trials

T=T~=0;

T=tau(i).*T;

FracTake=SIS_frac_sim(T,g,No_sims);

Frac_take(:,i)=[tau(i); FracTake];

end

Frac_take

toc

profile off

function FracTake=SIS_frac_sim(T,g,No_sims)

Fail_Vec=[]; I0=randsample(length(T(:,1)),1);

for i=1:No_sims

Fail=simulate(T,I0,g);

Fail_Vec=[Fail_Vec Fail];

FracTake=(length(Fail_Vec)-sum(Fail_Vec))/length(Fail_Vec);

end

function Fail=simulate(T,I0,g)

Fail=[]; No=length(T(:,1));

I_vec=zeros(1,No);I_vec(I0)=1;

S_vec=ones(1,No);S_vec(I0)=0;

S_tot=sum(S_vec);I_tot=sum(I_vec);

M=I_vec*T; P_vec=S_vec.*M; P=max(sum(P_vec),eps);

current_time=0; infection_time=current_time+exprnd(1/P);

R=max(g*sum(I_vec),eps); recovery_time=current_time+exprnd(1/R);

event=0; t_max=1000; I_max=20;

while I_tot<I_max

event=event+1;

if I_tot>0

%find next event

[time_to_next_event,event_type]=min([infection_time,recovery_time,t_max+1]);

current_time=time_to_next_event;

if event_type==1 %infection

A=P_vec/P;

edges=[0,cumsum(A)];

[F,farm_index]=histc(rand,edges);

S_vec(farm_index)=0; I_vec(farm_index)=1;

M=M+T(farm_index,:);

P_vec=max(S_vec.*M,0); P=max(sum(P_vec),eps);

infection_time=current_time+exprnd(1/P);

R=max(g*sum(I_vec),eps);

recovery_time=current_time+exprnd(1/R);

S_tot=S_tot-1; I_tot=I_tot+1;

elseif event_type==2 %recovery

A=I_vec/sum(I_vec);

edges=[0,cumsum(A)];

[F,farm_index]=histc(rand,edges);

I_vec(farm_index)=0; S_vec(farm_index)=1;

M=M-T(farm_index,:);

P_vec=max(S_vec.*M,0); P=max(sum(P_vec),eps);

infection_time=current_time+exprnd(1/P);

R=max(g*sum(I_vec),eps);

recovery_time=current_time+exprnd(1/R);

I_tot=I_tot-1; S_tot=S_tot+1;

end

else

Fail=1;

return

end

end

if I_tot==I_max; Fail=0; end

Reed-forst epidemic

n=66

p=.066

sick = 1;

SIR = zeros(1,n+1);

x = 0:n;

tau = 0;

healthy=n+1; %number of healthy people

while healthy > 1 && sick>0 %run while there are sick & healthy people

for Sz = 0:x(healthy); %potential number of people who don't get sick

SIR(1, Sz+1) = [nchoosek(x(healthy), Sz) * (1-p)^(sick*Sz) * (1-(1-p)^sick)^(x(healthy)-Sz)] ;

end

SIR(1, x(healthy)+2:end) = zeros(1,length(SIR(1, x(healthy)+2:end))); %delete old probabilities

cdf = cumsum(SIR(1,:));

u = rand;

healthy = 1;

while cdf(healthy) <= u %choose an new number of healthy people according to trans.prob

healthy = healthy+1;

end

x(healthy) %new number of healthy people in population +1

sick = size(find(SIR),2) - x(healthy) %new number of sick people

tau = tau + 1;

end

Related Questions

Navigate

• Browse (All)
• Browse S
• Subjects

---

---

Integrity-first tutoring: explanations and feedback only — we do not complete graded work. Learn more.