clear, close

p = 0.5;
num_trials = 100;
num_samples = 4000;
mean_vec = NaN(1,num_samples);

for k = 1:num_samples
Y_vec = NaN(1,num_trials);

for m = 1:num_trials

Y = 0;
c = 0;

X = NaN(1,7);

for n = 1:7
   r = rand;
   
   if r<p
       X(n) = 1;
   else
       X(n) = 0;
   end
   
   if X(n)==1
       c = c+1;
   else
       c = 0;
   end
   
   if c>=3
      Y = 1; 
   end
end

Y_vec(m) = Y;

end

rainy_weeks = sum(Y_vec);
mean_vec(k) = rainy_weeks;

end
avg = mean(mean_vec);

figure, hold on
histogram(mean_vec,30)
yL = get(gca,'YLim');
line([avg avg],yL,'Color','r')
xlabel('Number of rainy weeks out of 100')
ylabel('Count')

fprintf('Mean of sample distribution: %.2f \n',avg)