%% Problem Sheet 4
% Methods of Monte Carlo Simulation

clear all; close all; clc;

%% Exercise 2 (conditional expectation)

lambda = 5;
p = 0.2;
k = 2;

N = 10^4;

rejections = 0; % (not needed)
sample = zeros(N,1);
for i=1:N
    M = poissrnd(lambda);
    K = sum(rand(M,1) <= p);
    while ~(K == k)
        M = poissrnd(lambda);
        K = sum(rand(M,1) <= p);
        rejections = rejections + 1;
    end
    sample(i) = M;
end

fprintf('The estimated conditional expectation of M given K=%d is %f.\n', k, mean(sample));
fprintf('The estimated standard deviation of our estimator is %f.\n', std(sample)/sqrt(N));

% true value: (1-p)*lambda + k = 6  (see theory)

% not needed
fprintf('The acceptance rate is %f.\n', N / (N+rejections));


%% Exercise 5 (casino evening)

close all; clear all;

% a)

N = 10^5;

p = 0.1; % parameter of the geometric distribution we want to sample from
         % (we want the first failure, so p has to be the failure prob.)
lambda = -log(1-p); % parameter of the corresponding exponential distribution

X = zeros(N,1);
for i = 1:N
    Z = -(1/lambda)*log(rand());
    X(i) = floor(Z); % (or ceil(Z)-1, we don't want to count the beer we paid)
end

fprintf('The expected number of beers you can drink for free is %f\n', mean(X));


% b)

p = 0.5; % success probability (with which the bartender gets money)
lambda = -log(1-p); 
lambda2 = 30; % parameter of the Poisson distribution

Y = zeros(N,1);
for i = 1:N
    % Poisson number of people, K
    s = 0; K = 0;
    Z = - (1/lambda2)*log(rand());
    s = s + Z;
    while s <= 1
        K = K + 1;
        Z = - (1/lambda2)*log(rand());
        s = s + Z;
    end
    % Given K, generate the appropriate binomial
    m = 0; s = 0;
    while m <= K
        Z = -(1/lambda)*log(rand());
        V = ceil(Z); % this time we need the version counting the success
        m = m + V;
        s = s + 1;
    end
    Y(i) = s-1; % minus 1 because we counted one too many
end

figure(1);
counts = hist(Y, 0:max(Y));
bar(0:max(Y), counts/sum(counts), 1);
title('Distribution of the number of beers that are paid');
ylabel('relative frequency');


%% Exercise 7 (scandal)

close all; clear all;

N = 10^4;

% parameters of the negative binomial distribution
r = 6;
p = 0.8;
lambda = -log(1-p); % parameters of the suitable geometric distribution

X = zeros(N,1);
for i = 1:N
    n = 3; % initial number of witnesses
    for j = 1:5 % loop for days
        Y = zeros(n,1);
        for k = 1:n % loop for persons knowing about the scandal
            % generate a negative binomial like in Ex. 6
            Z = floor(-(1/lambda)*log(rand(r,1)));
            Y(k) = sum(Z);
        end
        % add new persons who know about the scandal
        % (because they help spreading the news the next day)
        n = n + sum(Y);
    end
    % end of one simulation run -> remember final number of people
    X(i) = n;
end

fprintf('On average, %f people know about the scandal after 5 days.\n', mean(X));