%% 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));