clear all
clf
%a) 
N=6;
I1=1;
s1=N-I1;
p=0.27;
nfamilies=100;
q=(1-p)^I1;
for(k=0:s1)
dist_s2(k+1)=nchoosek(s1,k).*q.^k.*(1-q).^(s1-k);
end
average_s2=s1*q
k=[0:s1];
figure(1)
bar(k,dist_s2)

%b) Version 1
%for each family
for(i=1:nfamilies)
    In=I1;
    q=(1-p)^I1;
    s(i)=s1;
%the disease can take a maximum of s1 days (until everybody either gets
%sick or not), go through every day
    for(time=1:s1)
    countS=0;
%    for each person who is still healthy
    for(l=1:s(i))
        %toss a coin, if still stays healthy, count the person
        r=rand(1);
        if(r<q)
            countS=countS+1;
        end
    end
    %final count
    In=s(i)-countS;
    s(i)=countS;
    q=(1-p)^In;
end
end
figure(2)
clf
plot(s,'o')

%b) version 2

for(j=1:nfamilies)
    S(j)=s1;
    I=I1;
    for(time=1:s1)
        q=(1-p).^I;
        for(k=0:S(j))
          dist(k+1)=nchoosek(S(j),k).*q.^k.*(1-q).^(S(j)-k);
        end
        cumdist=cumsum(dist);
        r=rand(1);
        nhealthy=min(find(cumdist>r))-1;
        I=S(j)-nhealthy;
        S(j)=nhealthy;
    end
end
hold on 
plot(S,'or')
          
figure(3)
clf
for(k=0:s1)
    bar(k,length(find(s==k)))
    hold on
end