
%%%%% SIR stochastic simulation

%% use "next event" simulation technique

sizes=zeros(1000,1);
for j=1:1000   % do 1000 realizations of the model

t=0;
n=1000;  
i=1;   % one initial infective

s=n-i;
tc=i;

beta=100.0;
gamma=50.0;

beta_over_n=beta/n;

while (i>0)
    t(1)=exprnd( 1/(beta_over_n*s*i) ); % when would next infection occur?
    t(2)=exprnd( 1/(gamma*i)  ); % when would next recovery occur?
    
     
% which happened first: infection or recovery?

    if (t(1)<t(2))    
           s=s-1;
           i=i+1;
           tc=tc+1;
           t=t+t(1);
       else       
           i=i-1;
           t=t+t(2);                  
    end
    
    
end
j,tc  % print out realization number and outcome
      % (mainly so that we can see our progress)
sizes(j)=tc;
end
hist(sizes,100)  
% the "100" means use 100 bars in the histogram


    