%% Problem Sheet 4 % Stochastic Simulation close all; clear all; clc; %% Exercise 3 (acceptance-rejection) n = 10; % parameters of the binomial distribution p = 0.2; N = 10^4; % sample size % a) (uniform distribution) f = binopdf(0:10, n, p); % binomial pdf g = ones(1, 11)/11; % uniform pdf C = max(f./g); % optimal scaling constant Y1 = []; while length(Y1) < N X = floor(11*rand()); % proposal U = rand(); if U <= f(X+1)/(C*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 a)'); figure(2); bar(0:10, f, 1); title('Binomial distribution (n=10, p=0.2)'); % b) p2 = 0.3; % parameter for geometric distribution lambda = -log(1-p2); % corresponding parameter for exponential distribution f = binopdf(0:10, n, p); % binomial pdf g = geopdf(0:10, p2); % geometric pdf (where needed) C = max(f./g); % optimal scaling constant Y2 = []; while length(Y2) < N X = floor(-log(rand())/lambda); U = rand(); if (X <= 10) && (U <= f(X+1)/(C*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 b)'); %% Exercise 4 % a) close all; clear all; n = 5; % (3, 5, 10) N = 10^3; success = zeros(N,1); rejections = 0; for i=1:N X = ceil(n*rand(1,n)); while max(hist(X, 1:n)) > 1 X = ceil(n*rand(1,n)); % count rejections rejections = rejections+1; end success(i) = (abs(find(X==1) - find(X==n)) == 1); end fprintf('The estimated acceptance probability is %f\n', N/(rejections+N)); fprintf('The probability of the cars standing next to each other is %f\n', mean(success)); % theoretical acceptance probability factorial(n)/(n^n) % theoretical probability for cars standing next to each other 2*(n-1)*factorial(n-2)/factorial(n) % b) close all; clear all; N = 10^3; n = 3; % (3, 5, 10) p = 0.1; success = zeros(N,1); rejections = 0; for i=1:N % determine which employees are ill ill = (rand(1,n) <= p); % acceptance-rejection for the positions X = floor((n+1)*rand(1,n)); % check validity h = hist(X,0:n); accept = (h(1) == sum(ill)) && (max(h(2:end)) <= 1); for j = 1:length(ill) if ill(j) && (h(j+1) > 0) accept = false; break; end end % if not, repeat simulation and check while ~accept X = floor((n+1)*rand(1,n)); h = hist(X,0:n); accept = (h(1) == sum(ill)) && (max(h(2:end)) <= 1); for j = 1:length(ill) if ill(j) && (h(j+1) > 0) accept = false; break; end end % count rejections rejections = rejections+1; end if (h(2)>0) && (h(n+1) > 0) success(i) = (abs(find(X==1) - find(X==n)) == 1); end end fprintf('The estimated acceptance probability is %f\n', N/(rejections+N)); fprintf('The probability of the cars standing next to each other is %f\n', mean(success)); %% Exercise 10 close all; clear all; % a) % State space: {0,1}^(5x5) % Proposal distribution: q_X,Y = 1/(n^2) if X and Y differ in exactly % 1 element (zero otherwise) % Acceptance probability: 1 if the configuration is legal, zero otherwise % b) N = 10^5; % sample size n = 5; % dimension X = zeros(n,n); % initial state for i=1:N % generate a proposal index = ceil(n*rand(1,2)); Y = X; Y(index(1),index(2)) = 1-Y(index(1),index(2)); % check if it is legal legal = true; for j=1:n for k=1:(n-1) if Y(j,k)==1 && Y(j,k+1)==1 legal = false; break; end if Y(k,j)==1 && Y(k+1,j)==1 legal = false; break; end end if ~legal break; end end % if yes, accept; otherwise reject and stay in X. if legal X = Y; end end X