%% Problem Sheet 9 % Stochastic Simulation close all; clear all; clc; %% Exercise 3 % initialization N = 10^4; X = zeros(N+1, 1); X(1) = 0; f = @(x) (exp(-(x-2).^2/2) + exp(-(x+2).^2/2))/(2*sqrt(2*pi)); for i = 1:N Y = rand() * f(X(i)); % shift function f such that g is zero if f(x)=Y g = @(x) f(x) - Y; % compute start/end points of interval(s) a = fzero(g, 1.5); b = fzero(g, 2.5); if abs(b-a)<10^-6 % found the same point twice -> one large interval X(i+1) = 2*a*rand() - a; else % two smaller intervals x = rand(); if x > 0.5 X(i+1) = 2*(x-0.5)*(b-a) + a; else X(i+1) = 2*(x-0.5)*(b-a) - a; end end end % plot figure(1); h = hist(X, -10:0.5:10); bar(-10:0.5:10, h/(sum(h)*0.5), 1); hold on; x = -10:0.1:10; plot(x, f(x), 'green'); hold off; axis([-10, 10, 0, 0.3]) %% Exercise 4 clear all; f = @(x) x.^2 .* exp(-x/2) .* (x >= 0); N = 10^5; X = zeros(N,1); X(1) = 2; for i = 2:N % proposal: current value + std normal Y = X(i-1) + randn(); % acceptance probability alpha = f(Y)/f(X(i-1)); if rand() < alpha X(i) = Y; else X(i) = X(i-1); end end % plot figure(2); h = hist(X, 0:50); bar(0:50, h/sum(h), 1); hold on; x = 0:0.1:50; C_est = 1/(sum(f(x))*0.1); plot(x, C_est*f(x), 'green'); hold off; axis([-10, 40, 0, 0.2]);