% This program runs the epidemic simulator lots of times and
% compares to some analytic results.

N=100 %number of individuals
T=2000 %time steps of simulation
R=100 %Number of replicates. Increasing R gives a more accurate outcome of the model

vq=0;
p=0.01;


range=[0.001:0.001:0.2];
for q=range   
    vq=vq+1
    for runs=1:R
        
        %run the simulation
        x=runepidemic(p,q,N,T);

        %Get the mean population over the last 500 time steps
        fx(vq,runs)=mean(x(T-500:T));
        %Get the population on the last time step
        distx(vq,runs)=x(T);

    end

end

figure(3)
plot(range(1:124),fx(1:124,:),'yx')
hold on
plot(range(1:124),N*(1-p./range(1:124)),'k')
xlabel('q')
ylabel('Equilibrium infected population')
axis([0 0.13 -1 100])

figure(4)
plot(range(1:124),sum(fx(1:124,1:R)'==0)/R,'x')
hold on
plot(range(1:124),p*(1-range(1:124))./(range(1:124).*(1-p)),'k')
xlabel('q')
ylabel('Proportion of cases where the disease spread')
axis([0 0.13 0 1.1])


hrange=[0:1:100];
histx=zeros(124,length(hrange));
for i=1:124
    histx(i,:)=hist(dx(i,:),hrange);
end

figure(5)
imagesc(range(1:124),hrange,histx'/R,[0 1])
xlabel('q')
ylabel('Innfected individuals in population')
colormap

