%% Problem Sheet 3 % Methods of Monte Carlo Simulation clear all; close all; clc; %% Exercise 1 (discrete acceptance-rejection) % a) and b) n = 10; % parameters of the binomial distribution p = 0.2; N = 10^4; % sample size % (i) (uniform distribution) f = binopdf(0:10, n, p); % binomial pdf g = ones(1, 11)/11; % uniform pdf C_i = max(f./g); % optimal scaling constant Y1 = []; while length(Y1) < N X = floor(11*rand()); % proposal U = rand(); if U <= f(X+1)/(C_i*g(X+1)) % condition for accepting Y1 = [Y1, X]; end end figure(1); counts = hist(Y1, 0:10); bar(0:10, counts/sum(counts), 1); title('Histogram of values simulated in (i)'); figure(2); bar(0:10, f, 1); title('Binomial distribution (n=10, p=0.2)'); % (ii) (geometric distribution) q = 0.3; % parameter for geometric distribution lambda = -log(1-q); % corresponding parameter for exponential distribution f = binopdf(0:10, n, p); % binomial pdf g = geopdf(0:10, q); % geometric pdf (where needed) C_ii = max(f./g); % optimal scaling constant Y2 = []; while length(Y2) < N % generate exponential using inverse transform X = floor(-log(rand())/lambda); U = rand(); if (X <= 10) && (U <= f(X+1)/(C_ii*g(X+1))) Y2 = [Y2, X]; end end figure(3); counts = hist(Y2, 0:10); bar(0:10, counts/sum(counts), 1); title('Histogram of values simulated in (ii)'); % c) fprintf('The acceptance probabilities are %f for (i) and %f for (ii).\n', 1/C_i, 1/C_ii); fprintf('=> The geometric proposals should be preferred.\n'); %% Exercise 3 (sampling from a shape) close all; clear all; R = 5; r = 2; N = 10^4; rejections = 0; X = zeros(N,3); for i = 1:N Y(1:2) = rand(1,2)*2*(r+R) - (r+R); Y(3) = rand()*2*r - r; while (R - sqrt(Y(1)^2 + Y(2)^2))^2 + Y(3)^2 > r^2 Y(1:2) = rand(1,2)*2*(r+R) - (r+R); Y(3) = rand()*2*r - r; rejections = rejections+1; end X(i,:) = Y; end figure(1); plot3(X(:,1), X(:,2), X(:,3), 'o'); axis([-10 10 -10 10 -10 10]); title('Points uniformly distributed in a torus'); fprintf('The estimated acceptance probability is %f\n', N/(N+rejections)); %% Exercise 6 (Box-Muller) close all; clear all; % a) % see function box_muller N = 10^5; X = box_muller(N); % histogram (normalized, so the area of bars is 1) figure(1); w = 0.5; % width of bins h = hist(X, -20:w:20); bar(-20:w:20, h/sum(h)/w, 1); axis([-7, 7, 0, 0.6]); ylab('relative frequency ()'); % add density hold on; x = linspace(-10, 10, 100); plot(x, (1/sqrt(2*pi))*exp(-x.^2/2), 'g'); hold off; % b) mu = 20; sigma_sq = 100; X = sqrt(sigma_sq)*box_muller(N) + mu; figure(2); hist(X, -40:5:80); axis([-30, 70, 0, 3e4]); ylab('absolute frequency'); title(sprintf('Histogram of N(%d,%d)-distributed random numbers', mu, sigma_sq));