CMPE58N - Monte Carlo MethodsExamples covered on 14.05Consider the following dynamical sytem. Design two SMC methods, first one by changing the observation model and defining a new likelihood function, and the second one using the exact model and a smart proposal distribution Particle Filtering assignment
% function [] = cmpe58n_rotor_pf2() % CMPE58N_ROTOR_PF2 Temporary Script created 14-May-2014 % % [] = CMPE58N_ROTOR_PF2() % % Inputs : % : % % Outputs : % : % % Usage Example : [] = cmpe58n_rotor_pf2(); % % % Note : % See also % Uses : % Change History : % Date Time Prog Note % 14-May-2014 11:28 AM ATC Created under MATLAB 8.0.0.783 (R2012b) % ATC = Ali Taylan Cemgil, % Department of Computer Engineering, Bogazici University % e-mail : taylan.cemgil@boun.edu.tr rotmat = @(th)[cos(th) -sin(th);sin(th) cos(th)]; theta = 2*pi/100; A = blkdiag(rotmat(theta), rotmat(2*theta)); C = [1 0 0 0;0 0 0 1]; x = [1 0 1 0]'; T = 100; y = zeros(2, T); % Sensor locations L = [1 1;-1 -1]'; % Number of sensors M = size(L,2); z = zeros(M, T); Q = 0.001; R = 0.001; for t=1:T, y(:,t) = C*x + sqrt(R)*randn(2,1); x = A*x + sqrt(Q)*randn(4,1); for m=1:M, z(m, t) = norm(y(:,t)-L(:,m)); end end subplot(211) plot(y(1,:), y(2,:), '-o'); axis equal subplot(212) plot(z') Examples covered on 12.03Importance Sampling Illustration 2
% CMPE58N_IS02 Estimate the tail probability of a Gaussian % Change History : % Date Time Prog Note % 12-Mar-2014 11:17 AM ATC Created under MATLAB 7.10.0.499 (R2010a) % ATC = Ali Taylan Cemgil, % Department of Computer Engineering, Bogazici University % e-mail : taylan.cemgil@boun.edu.tr %% Estimate the probability Pr{x>L} where x ~ N(x; 0, 1) L = 8; qfunc(L) %% vanilla MC estimate % Number of samples N = 100000; x = randn(N, 1); % test function phi = @(x, L)double(x > L); p_est = 1/N*sum(phi(x, L)) %% IS % log of the target density N(x; 0, 1) lp = @(x)-0.5*log(2*pi)-0.5*x.^2; % log of the proposal density N(x; mu, sig2) mu = L; sig2 = 2; lq = @(x)-0.5*log(2*pi*sig2)-0.5*(x-mu).^2./sig2; % Generate samples from the proposal x = sqrt(sig2)*randn(N, 1) + mu; % Calculate weights lw = lp(x) - lq(x); % Stable calculation of log-weights lw2 = lw - max(lw); w = exp(lw2)/sum(exp(lw2)); % IS Estimate p_est2 = sum(w.*phi(x, L)) Importance Sampling Illustration 1
% CMPE58N_IS01 Estimate the variance of a Beta RV % % Change History : % Date Time Prog Note % 12-Mar-2014 10:18 AM ATC Created under MATLAB 7.10.0.499 (R2010a) % ATC = Ali Taylan Cemgil, % Department of Computer Engineering, Bogazici University % e-mail : taylan.cemgil@boun.edu.tr %% Parameters of the beta distribution a = 40; b = 80; v_true = a*b/((a+b)^2*(a+b+1)) %% MC: % Number of samples N = 1000; % Number of experiment E = 100; % Variance estimate at each experiment V_EST1 = []; for e = 1:E, x = betarnd(a,b, N, 1); v_est1 = 1/N*sum(x.^2) - (1/N*sum(x)).^2; V_EST1 = [V_EST1 v_est1]; end; subplot(211) hist(V_EST1, 20); xlim = [0.0 0.005]; set(gca, 'xlim', xlim); %% IS % target density upto a normalization phi = @(x, a, b)x.^(a-1).*(1-x).^(b-1); % Proposal is the uniform density on [0 1] % Normalization constant, we don't need this Z = @(a, b)gamma(a+b)/(gamma(a)*gamma(b)); V_EST2 = []; for e = 1:E, x = rand(N, 1); W = phi(x, a, b)/1; w = W./sum(W); v_est2 = sum(w.*x.^2) - (sum(w.*x)).^2; V_EST2 = [V_EST2 v_est2]; end subplot(212) hist(V_EST2, 20); set(gca, 'xlim', xlim); Examples from 2012Examples covered on 08.05Particle Filter for the Stochastic Volatility Model
% CMPE58N_SV_PF PF Demo % % [] = CMPE58N_SV_PF() % % Inputs : % : % % Outputs : % : % % Usage Example : [] = cmpe58n_sv_pf(); % % % Note : % See also % Uses : % Change History : % Date Time Prog Note % 08-May-2013 12:15 PM ATC Created under MATLAB 7.10.0.499 (R2010a) % ATC = Ali Taylan Cemgil, % Department of Computer Engineering, Bogazici University % e-mail : taylan.cemgil@boun.edu.tr % cmpe58n_sv_model N = 1000; % number of particles T = model.T; y = model.y; sigma = model.sigma; npdf = @(y,l)exp(-0.5*y.^2./exp(l))./sqrt(2*pi*exp(l)); l0 = zeros(1, N); W = ones(1,N)/N; l_est = zeros(1, T); for t=1:T, % Propose if t==1, l_new = l0 + sigma*randn(size(l0)); else l_new = l + sigma*randn(size(l)); end W = W.*npdf(y(t), l_new); w = W./sum(W); W = w; l_est(t) = w*l_new'; ESS = 1./sum(w.^2); ESS % if (ESS < N/2) % Resample if 0, [idx] = cmpe58n_resample(w); l = l_new(idx); W = ones(1,N)/N; else l = l_new; end; end; plot(y) hold on plot(3*exp(model.l_true/2), 'r:'); plot(3*exp(l_est/2), 'k'); hold off Stochastic Volatility Model
% CMPE58N_SV_MODEL Temporary Script created 08-May-2013 % % [] = CMPE58N_SV_MODEL() % % Inputs : % : % % Outputs : % : % % Usage Example : [] = cmpe58n_sv_model(); % % % Note : % See also % Uses : % Change History : % Date Time Prog Note % 08-May-2013 11:57 AM ATC Created under MATLAB 7.10.0.499 (R2010a) % ATC = Ali Taylan Cemgil, % Department of Computer Engineering, Bogazici University % e-mail : taylan.cemgil@boun.edu.tr sigma = 0.06; l0 = 0; T = 5000; l = zeros(1, T); y = zeros(1, T); for t=1:T, if (t>1) l(t) = l(t-1) + sigma*randn; else l(1) = l0 + sigma*randn; end y(t) = exp(l(t)/2)*randn; end plot(filter([ones(1, 10)], 1, abs(y))) hold on plot(3*exp(l/2), 'r'); hold off model = struct('sigma', sigma, 'y', y, 'T', T, 'l_true', l) Systematic Resampling
function [idx] = cmpe58n_resample(W) % CMPE58N_RESAMPLE Systematic Resampling of a normalized weight vector % % [idx] = CMPE58N_RESAMPLE(W) % % Inputs : % : % % Outputs : % : % % Usage Example : [] = cmpe58n_resample(); % % % Note : % See also % Uses : % Change History : % Date Time Prog Note % 08-May-2013 11:54 AM ATC Created under MATLAB 7.10.0.499 (R2010a) % ATC = Ali Taylan Cemgil, % Department of Computer Engineering, Bogazici University % e-mail : taylan.cemgil@boun.edu.tr N = length(W); u = rand/N + ((1:N)-1)/N; [~, idx] = histc(u, [0 cumsum(W)]); |