%Stochastic Reed-Frost Model

ngen=20; %number of generations to simulate
Sinit=2000; %initial susceptible population
Iinit=2; %initial infected population
Rinit=0; %initial recovered/immune population

q=0.999;

nsims=10; %number of simulations to run

for i=1:nsims  
    S=zeros(ngen,1);
    I=zeros(ngen,1);
    R=zeros(ngen,1);
    S(1)=Sinit;
    I(1)=Iinit;
    R(1)=Rinit;
    for j=1:ngen
        I(j+1)=binornd(S(j),1-q^I(j));
        S(j+1)=S(j)-I(j+1);
        R(j+1)=R(j)+I(j);
    end;
    plot(0:ngen,I,'*-')
    ylabel('Number of Infectives');
    xlabel('Generation');
    hold on
end; %simulation loop

    



