%% Problem Sheet 3 % Stochastic Simulation close all; clear all; clc; %% Exercise 5 N = 30; % persons n = [100, 1000, 10000]; % steps m = max(n); P = zeros(N+1, N+1); for i = 1:N P(i, i+1) = (N-(i-1))/N; end for i = 2:(N+1); P(i, i-1) = (i-1)/N; end % simulation of MC X = zeros(m, 1); X(1) = 0; for j = 1:m U = rand; % states and indices are shifted -> +-1 X(j+1) = find(cumsum(P(X(j)+1,:)) > U, 1) - 1; end % histograms for i = 1:3 figure(i); hist(X(1:n(i)), 0:30); end % %%% not asked in exercise % % % nice to know: no limiting distribution % % the Markov chain is periodic! % % m_even = 1:2:m; % m_odd = 2:2:m; % figure(4); % hist(X(m_even), 0:30); % figure(5); % hist(X(m_odd), 0:30); % %%% %% Exercise 6 close all; clear all; % transition matrix and invariant distribution of Y_n P = [0, 3/4, 1/4; 1/2, 0, 1/2; 1/3, 1/3, 1/3]; mu = [12/41, 14/41, 15/41]; % stationary initial distribution % transition matrix of time-reversal P_hat = zeros(3,3); for i=1:3 for j=1:3 P_hat(i,j) = P(j,i)*mu(j)/mu(i); end end % a) n = 50; X = zeros(n+1,1); % simulate X_n, i.e., X_0 of the backwards chain X(1) = find(rand() < cumsum(mu), 1 ); for i=1:n % simulate each step using P_hat X(i+1) = find(rand() < cumsum(P_hat(X(i),:)), 1); end figure(1); clf for i=1:n line([i-1, i], [X(i), X(i)]); end axis([0, n, -1, 5]); % b) n = 10^5; X = zeros(n+1,1); P_est = zeros(3,3); X(1) = find(rand() < cumsum(mu), 1); for i=1:n X(i+1) = find(rand() < cumsum(P_hat(X(i),:)), 1); % count (backwards) transitions P_est(X(i+1),X(i)) = P_est(X(i+1),X(i))+1; end % normalize P_est for i=1:3 P_est(i,:) = P_est(i,:)/sum(P_est(i,:)); end P_est P