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