%% Problem Sheet 1 % Methods of Monte Carlo Simulation %% Exercise 1 (surface area estimation) close all; clear all; clc; % circle centers and radius center1 = [0, 0.5]; center2 = [1, 0.5]; r = 1/sqrt(2); % sample sizes N = [10^2, 10^3, 10^5]; for n = N % initialize vector for counting how many points fell inside insideA = zeros(n,1); for i = 1:n % draw random point in (0,1)^2 x = rand(1,2); % check if it is in the gray region insideA(i) = (sum((center1-x).^2) <= r^2) && (sum((center2-x).^2) <= r^2); end fprintf('The estimated area using %d points is: %f\n', n, mean(insideA)); end fprintf('The true surface area is: %f\n', pi/4 - 1/2); %% Exercise 2 (RANDU) close all; clear all; % sample size N = 10^4; X = zeros(N,1); S = zeros(N+1,1); % seed S(1) = 15; % RANDU a = 2^16 + 3; c = 0; m = 2^31; for i = 1:N S(i+1) = mod(a*S(i) + c, m); X(i) = S(i+1) / m; end % 1-dim plot figure(1) plot(X); % 2-dim plot figure(2) scatter(X(1:N-1), X(2:N)); % 3-dim plot figure(3) scatter3(X(1:N-2), X(2:N-1), X(3:N)); %% Exercise 3 (Monte Carlo Integration) close all; clear all; % sample size N = 10^4; M = 3 * N; U = zeros(M, 1); V = zeros(M+3, 1); W = zeros(M+3, 1); % parameters of generator m1 = 2^32 - 209; m2 = 2^32 - 22853; a11 = 1403580; a12 = 810728; a21 = 527612; a22 = 1370589; % seed V(1) = 1359; V(2) = 4976; V(3) = 1279; W(1) = 5879; W(2) = 8265; W(3) = 2269; for i = 1:M V(i+3) = mod(a11 * V(i+1) - a12 * V(i), m1); W(i+3) = mod(a21 * W(i+2) - a22 * W(i), m2); if V(i+3) <= W(i+3) U(i) = (V(i+3) - W(i+3) + m1) / (m1 + 1); else U(i) = (V(i+3) - W(i+3)) / (m1 + 1); end end % split into 3 vectors X1 = U(1:3:end); X2 = U(2:3:end); X3 = U(3:3:end); % plug into function Y = cos(X3.*X2) .* sin(cos(X3.*X1)); int_est = mean(Y); fprintf('The estimated value of the integral is: %f\n', int_est);