%% Problem Sheet 6 % Methods of Monte Carlo Simulation clear all; close all; clc; %% Exercise 2 (estimation of pi) N = 10^4; % a) U_1 = rand(N,1); X_a = 4*sqrt(1 - U_1.^2); fprintf('The estimate of pi using conditional MC is %f.\n', mean(X_a)); % b) U = rand(N/2, 2); Y = 4*(U(:,1).^2 + U(:,2).^2 <= 1); Y_star = 4*((1-U(:,1)).^2 + (1-U(:,2)).^2 <= 1); X_b = [Y; Y_star]; fprintf('The estimate of pi using antithetic variables is %f.\n', mean(X_b)); % Justification of the estimator: % % If (U1, U2) is in the quarter of the circle, (1-U1, 1-U2) is likely to be % not. The other way round, if (U1, U2) is not in the circle, (1-U1, 1-U2) % surely is. % % => The two indicators have to be negatively correlated, the constant does % not change this. % c) fprintf('The estimated variances of the estimators are %e for a) and %e for b).\n\n', ... var(X_a)/N, var(X_b)/N); % => Conditional MC performs better in this case. %% Exercise 4 (random walk) close all; clear all; N = 10^4; % a) (crude Monte Carlo) % generate Z_i U = rand(N, 10); Z = zeros(N,10); Z(U < 0.5) = 1; Z(U >= 0.5) = -1; % check if sums are > 1 Y = (sum(Z, 2)>1); fprintf('Estimated probability using standard MC: %f\n', mean(Y)); % b) % generate Z_i for N/2 random walks U2 = rand(N/2, 10); Z2 = zeros(N/2, 10); Z2(U2 < 0.5) = 1; Z2(U2 >= 0.5) = -1; % check if X_10 (the sum of Z_1, ..., Z_10) is > 1 ... Y2 = (sum(Z2, 2) > 1); % ... for the original random walks Y2_star = (-sum(Z2, 2) > 1); % ... and for their reflections fprintf('Estimated probability using antithetic variables: %f\n\n', mean([Y2; Y2_star])); % c) (Comparison) var_cmc = var(Y)/N; var_av = var(Y2 + Y2_star)/(2*N); fprintf('Variance of the estimator using standard MC: %d\n', var_cmc); fprintf('Variance of the estimator using antithetic variables: %d\n', var_av); fprintf('Quotient of the two: %f\n\n', var_av/var_cmc); %% Exercise 6 (stratified sampling) strata = [0, 1; 1, 2; 2, 6; 6, 10]; f = @(x) (9/50)*x.^3 -(7/2)*x.^2 + 20*x -22; % a) N = 10^3; U = rand(N,1)*10; p = zeros(size(strata,1), 1); sigma = zeros(size(strata,1), 1); for i=1:size(strata,1) X = U((U >= strata(i,1)) & (U < strata(i,2))); p(i) = length(X) / N; sigma(i) = std(10*f(X)); end q = sigma.*p / sum(sigma.*p); fprintf('The estimated proportions are: %f, %f, %f, %f\n', q); fprintf(['With a total sample size of 10^5, this corresponds to sample', ... ' sizes of: %d, %d, %d, %d\n\n'], round(q*10^5)); % b) N = 10^5; N_i = round(q .* N); Y_bar = zeros(size(strata,1),1); sigma_sq = zeros(size(strata,1),1); for i=1:size(strata,1) U = rand(N_i(i),1)*(strata(i,2)-strata(i,1)) + strata(i,1); Y_bar(i) = mean(10*f(U)); sigma_sq(i) = var(10*f(U))/N_i(i); end est_int = sum(p.*Y_bar); est_var = sum((p.^2).*sigma_sq); fprintf('---Stratified Sampling---\n'); fprintf('The estimated value of the integral is %f.\n', est_int); fprintf('The estimated variance is %f.\n\n', est_var); % c) U = rand(N,1)*10; est_int_2 = mean(10*f(U)); est_var_2 = var(10*f(U))/N; fprintf('---Standard Monte Carlo---\n'); fprintf('The estimated value of the integral is %f.\n', est_int_2); fprintf('The estimated variance is %f.\n\n', est_var_2);