%% Problem Sheet 5
% Methods of Monte Carlo Simulation

clear all; close all; clc;

% set initial seed
S = [sum(100*clock), sum(37*clock), sum(179*clock)];

%% Exercise 4

% very rough check of myrand
[x, S] = myrand(10^5, 1, S);
figure(1);
hist(x);
title('Histogram of values generated by myrand');

% rough check of mybetarnd
[y, S] = mybetarnd(10^5, 1, 2, 3, S);
h = 0.05;
counts = hist(y, 0:h:1);
figure(2);
bar(0:h:1, counts/(sum(counts)*h), 1);
hold on;
x = linspace(0, 1, 100);
plot(x, betapdf(x, 2, 3), 'g-');
hold off;
c = title('Histogram of values generated by mybetarnd for p=2 and q=3');
set(c, 'FontName', 'Times');

% c) 

N = 10^4; % sample size
kids = 2:6; % number of kids in each house
m = length(kids); % number of houses
p = 7; q = 2; % parameters of the beta distribution

kids_without_gift = zeros(N,1);
house_without_gift = zeros(N,1);

for i=1:N % loop for sample
    gifts = zeros(m,1);
    
    for j=1:m % loop for houses
        % draw gift probability
        [P, S] = mybetarnd(1, 1, p, q, S);
        % sample from appropriate binomial distribution
        [X, S] = myrand(kids(j), 1, S);
        gifts(j) = sum(X <= P);
    end
    
    kids_without_gift(i) = sum(kids) - sum(gifts);
    house_without_gift(i) = (min(gifts) == 0);
end

fprintf('The estimated expected number of kids without gift is %.2f.\n', ...
        mean(kids_without_gift));
fprintf('The estimated probability of a house without gifts is %f.\n', ...
        mean(house_without_gift));

fprintf('The estimated standard errors of the estimators are %f and %f.\n', ...
        std(kids_without_gift)/sqrt(N), ...
        std(house_without_gift)/sqrt(N));