Skip to content

Commit d1b646e

Browse files
author
compops
committed
added matlab skeletons
Former-commit-id: 47e42c8 [formerly 47e42c8 [formerly c4491a5]] Former-commit-id: 0686dd8c1dfdfa3902fe347729377a258b9ffd4c Former-commit-id: e89e1a9
1 parent fdb715b commit d1b646e

19 files changed

Lines changed: 1289 additions & 0 deletions

‎matlab-seminar/complete/earthquake_data.mat‎

494 Bytes
Binary file not shown.

‎matlab-seminar/complete/example1_lgss_state.m‎

Lines changed: 72 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,72 @@
1+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2+
%
3+
% Example of particle filtering for the linear Gaussian state space model
4+
%
5+
% Copyright (C) 2015 Johan Dahlin < johan.dahlin (at) liu.se >
6+
%
7+
% This program is free software; you can redistribute it and/or modify
8+
% it under the terms of the GNU General Public License as published by
9+
% the Free Software Foundation; either version 2 of the License, or
10+
% (at your option) any later version.
11+
%
12+
% This program is distributed in the hope that it will be useful,
13+
% but WITHOUT ANY WARRANTY; without even the implied warranty of
14+
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15+
% GNU General Public License for more details.
16+
%
17+
% You should have received a copy of the GNU General Public License along
18+
% with this program; if not, write to the Free Software Foundation, Inc.,
19+
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
20+
%
21+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22+
23+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24+
% Load data
25+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26+
load('lgss_data.mat');
27+
28+
% Set the parameters of the model (phi, sigmav, sigmae)
29+
par = [ 0.5, 1.0, 0.1 ];
30+
T = 250;
31+
32+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33+
% State estimation
34+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35+
36+
% Run the Kalman filter
37+
x0 = 0;
38+
P0 = 0.1;
39+
xhatf_kf = kf( y, par, T, x0, P0 );
40+
41+
% Run the particle filter
42+
% No. particles in the particle filter ( choose nPart ~ T )
43+
nPart = 100;
44+
[ xhatf_pf, ll_pf ] = sm_lgss( y, par, nPart, T, x0 );
45+
46+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
47+
% Plot the results
48+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
49+
50+
figure(1);
51+
% Plot the true state
52+
subplot(3,1,1);
53+
plot( 0:T, x );
54+
xlabel('time');
55+
ylabel('true state');
56+
57+
% Plot the state estimates
58+
subplot(3,1,2);
59+
plot( 0:T, xhatf_kf, 0:T, xhatf_pf );
60+
xlabel('time');
61+
ylabel('state estimate');
62+
legend('Kalman filter','Particle filter');
63+
64+
% Plot the squared error
65+
subplot(3,1,3);
66+
plot(0:T,log( (xhatf_kf - xhatf_pf).^2 ) );
67+
xlabel('time');
68+
ylabel('logarithm of squared error');
69+
70+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
71+
% End of file
72+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

‎matlab-seminar/complete/example2_earthquake_state.m‎

Lines changed: 61 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,61 @@
1+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2+
%
3+
% Example of particle filtering for the Earthquake model
4+
%
5+
% Copyright (C) 2015 Johan Dahlin < johan.dahlin (at) liu.se >
6+
%
7+
% This program is free software; you can redistribute it and/or modify
8+
% it under the terms of the GNU General Public License as published by
9+
% the Free Software Foundation; either version 2 of the License, or
10+
% (at your option) any later version.
11+
%
12+
% This program is distributed in the hope that it will be useful,
13+
% but WITHOUT ANY WARRANTY; without even the implied warranty of
14+
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15+
% GNU General Public License for more details.
16+
%
17+
% You should have received a copy of the GNU General Public License along
18+
% with this program; if not, write to the Free Software Foundation, Inc.,
19+
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
20+
%
21+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22+
23+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24+
% Load data
25+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
26+
load('earthquake_data.mat');
27+
28+
% Set the parameters of the model (mu, phi, beta)
29+
par = [ 0.88 0.15 16.58 ];
30+
T = 114;
31+
32+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
33+
% State estimation
34+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
35+
36+
% No. particles in the particle filter ( choose nPart ~ T )
37+
nPart = 100;
38+
39+
% Run the particle filter
40+
[ xhatf_pf, ll_pf ] = sm_earthquake( y, par, nPart, T );
41+
42+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
43+
% Plot the results
44+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
45+
46+
figure(2);
47+
% Plot the true state
48+
subplot(2,1,1);
49+
plot( 1:T, y );
50+
xlabel('time');
51+
ylabel('no. earthquakes');
52+
53+
% Plot the state estimates
54+
subplot(2,1,2);
55+
plot( 0:T, xhatf_pf );
56+
xlabel('time');
57+
ylabel('log-intensity');
58+
59+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
60+
% End of file
61+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

‎matlab-seminar/complete/example3_earthquake_parameters.m‎

Lines changed: 104 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,104 @@
1+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2+
%
3+
% Example of Particle Metropolis-Hastings (PMH) for the Earthquake model
4+
%
5+
% Copyright (C) 2015 Johan Dahlin < johan.dahlin (at) liu.se >
6+
%
7+
% This program is free software; you can redistribute it and/or modify
8+
% it under the terms of the GNU General Public License as published by
9+
% the Free Software Foundation; either version 2 of the License, or
10+
% (at your option) any later version.
11+
%
12+
% This program is distributed in the hope that it will be useful,
13+
% but WITHOUT ANY WARRANTY; without even the implied warranty of
14+
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15+
% GNU General Public License for more details.
16+
%
17+
% You should have received a copy of the GNU General Public License along
18+
% with this program; if not, write to the Free Software Foundation, Inc.,
19+
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
20+
%
21+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
22+
23+
24+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
25+
% Load data
26+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
27+
load('earthquake_data.mat');
28+
T = 114;
29+
30+
31+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
32+
% Parameter estimation using PMH
33+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34+
35+
% The inital guess of the parameter
36+
initPar = [ 0.5 0.5 15 ];
37+
38+
% No. particles in the particle filter ( choose nPart ~ T )
39+
nPart = 100;
40+
41+
% The length of the burn-in and the no. iterations of PMH
42+
% ( nBurnIn < nIter )
43+
nBurnIn = 250;
44+
nIter = 2000;
45+
46+
% The covariance matrix in the random walk proposal
47+
stepSize = diag( [0.07 0.03 2].^2 );
48+
stepSize = 0.8 * stepSize;
49+
50+
% Run the PMH algorithm
51+
[th, xh] = pmh_earthquake( y, initPar, nPart, T, nIter, stepSize );
52+
53+
% Compute posterior means
54+
thhat = mean( thhat, 1 )
55+
xhhat = mean( xh( nBurnIn:nIter, 2:(T+1) ), 1 )
56+
57+
58+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
59+
% Plot the results
60+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
61+
62+
figure(2);
63+
64+
% Plot the expected versus observed no. earthquakes
65+
subplot(4,2,[1 2]);
66+
plot( y )
67+
xlabel( 'time' );
68+
ylabel( 'no. earthquakes' );
69+
70+
hold on;
71+
plot( thhat(3) * exp(xhhat), 'r' );
72+
hold off;
73+
74+
% Plot the parameter posterior estimate
75+
% Plot the trace of the Markov chain after burn-in
76+
subplot(4,2,3);
77+
hist( th( nBurnIn:nIter, 1 ), floor( sqrt( nIter - nBurnIn ) ) );
78+
xlabel('phi');
79+
ylabel('posterior density estimate');
80+
81+
subplot(4,2,4);
82+
plot( nBurnIn:nIter, th( nBurnIn:nIter, 1 ) );
83+
xlabel('iteration'); ylabel('trace of phi');
84+
85+
subplot(4,2,5);
86+
hist( th( nBurnIn:nIter, 2 ), floor( sqrt( nIter - nBurnIn ) ) );
87+
xlabel('sigmav'); ylabel('posterior density estimate');
88+
89+
subplot(4,2,6);
90+
plot( nBurnIn:nIter, th( nBurnIn:nIter, 2 ) );
91+
xlabel('iteration'); ylabel('trace of sigmav');
92+
93+
subplot(4,2,7);
94+
hist( th( nBurnIn:nIter, 3 ), floor( sqrt( nIter - nBurnIn ) ) );
95+
xlabel('beta'); ylabel('posterior density estimate');
96+
97+
subplot(4,2,8);
98+
plot( nBurnIn:nIter, th( nBurnIn:nIter, 3 ) );
99+
xlabel('iteration'); ylabel('trace of beta');
100+
101+
102+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
103+
% End of file
104+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

‎matlab-seminar/complete/kf.m‎

Lines changed: 59 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,59 @@
1+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
2+
%
3+
% Example of Kalman filtering for the linear Gaussian state space model
4+
%
5+
% Bootstrap particle filter
6+
%
7+
% Copyright (C) 2015 Johan Dahlin < johan.dahlin (at) liu.se >
8+
%
9+
% This program is free software; you can redistribute it and/or modify
10+
% it under the terms of the GNU General Public License as published by
11+
% the Free Software Foundation; either version 2 of the License, or
12+
% (at your option) any later version.
13+
%
14+
% This program is distributed in the hope that it will be useful,
15+
% but WITHOUT ANY WARRANTY; without even the implied warranty of
16+
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17+
% GNU General Public License for more details.
18+
%
19+
% You should have received a copy of the GNU General Public License along
20+
% with this program; if not, write to the Free Software Foundation, Inc.,
21+
% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
22+
%
23+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
24+
25+
function xhatf = kf(y,par,T,x0,P0)
26+
27+
% Initalise variables
28+
xhatf = x0 * ones( T+1, 1 );
29+
xhatp = x0 * ones( T+1, 1 );
30+
Pp = P0;
31+
32+
% Set the parameters for the Kalman filter
33+
A = par(1);
34+
C = 1;
35+
Q = par(2)^2;
36+
R = par(3)^2;
37+
38+
% Run main loop
39+
for tt = 2:(T+1)
40+
41+
% Compute Kalman Gain
42+
S = C * Pp * C + R;
43+
K = Pp * C / S;
44+
45+
% Compute state estimate
46+
xhatf(tt) = xhatp(tt) + K * ( y(tt-1) - C * xhatp(tt) );
47+
xhatp(tt+1) = A * xhatf(tt);
48+
49+
% Update covariance
50+
Pf = Pp - K * S * K;
51+
Pp = A * Pf * A + Q;
52+
53+
end
54+
55+
end
56+
57+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
58+
% End of file
59+
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

‎matlab-seminar/complete/lgss_data.mat‎

4.08 KB
Binary file not shown.

0 commit comments

Comments
 (0)