CMPE58N - Monte Carlo Methods

Examples covered on 14.05

Consider 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.03

Importance 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 2012

Examples covered on 08.05

Particle 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)]);