%% Problem Sheet 10 % Stochastic Simulation close all; clear all; clc; %% Exercise 5 a) lambda = 5; t_max = 10; % 1. S = []; S_sum = 0; while S_sum < t_max % simulate inter-arrival times S_i = -log(rand())/lambda; S = [S S_i]; S_sum = S_sum + S_i; end % compute jump times T = cumsum(S); % for plot: add zero and truncate to t_max T = [0 T]; T(end) = t_max; N_t = 0:length(S)-1; figure(1); clf for i = 1:length(T)-1 line([T(i), T(i+1)], [N_t(i), N_t(i)]); end % 2. % make a grid for the time h = 0.05; t = 0:h:t_max; % simulate jump times (approximately) T = (rand(1, length(t)-1) <= lambda*h); % compute values of N_t on the grid N_t = [0 cumsum(T)]; figure(2); clf for i = 0:max(N_t) line([min(t(N_t==i)), max(t(N_t==i))], [i i]); end %% Exercise 5 b) N = 10^4; N_1_1 = zeros(N, 1); N_1_2 = zeros(N, 1); % make a grid for the time (only once) h = 0.05; t = 0:h:t_max; for i = 1:N % 1. S = []; S_sum = 0; while S_sum < t_max % simulate inter-arrival times S_i = -log(rand())/lambda; S = [S S_i]; S_sum = S_sum + S_i; end % compute jump times T = cumsum(S); N_1_1(i) = sum(T <= 1); % 2. % simulate jump times (approximately) J = (rand(1, length(t)-1) <= lambda*h); % compute values of N_t on the grid N_t_2 = [0 cumsum(J)]; N_1_2(i) = N_t_2(1/h + 1); end figure(1); h = hist(N_1_1, 0:20); bar(0:20, h/sum(h), 1); figure(2); h = hist(N_1_2, 0:20); bar(0:20, h/sum(h), 1); figure(3); bar(0:20, poisspdf(0:20, lambda), 1); %% Exercise 6 clear all; close all; f = @(x) 0.5*x.^2-10*exp(-(x.^2)/20).*cos(3*x); x = linspace(-10, 10, 500); % figure(1); % plot(x, f(x)); X_0 = 10; T_0 = 10; N = 10^3; sigma = 0.1; % [0.1, 0.5, 1] beta = 0.99; X = zeros(N+1,1); X(1) = X_0; T = T_0; for i=1:N Y = X(i) + sigma*randn(); alpha = exp( (f(X(i))-f(Y)) / T ); if rand() <= alpha X(i+1) = Y; else X(i+1) = X(i); end T = beta*T; end X(end) figure(2); [ax, p1, p2] = plotyy(x, f(x), X, 0:N); title('Simulated Annealing') xlabel('x'); ylabel(ax(1), 'f(x)'); ylabel(ax(2), 'iteration'); axis(ax(1), [-10 10 -20 60]); axis(ax(2), [-10 10 -100 N+300]) axis(ax(2), 'ij');