removed superfluous comments to improve readability. some minor refactorisations in matlab code.

Former-commit-id: 66e4794 [formerly 66e4794 [formerly ca86603]]
Former-commit-id: d32092f872e797c0b4dc71e7818052cef96d8ed9
Former-commit-id: 5b7efe5

Author: compops

.gitignore

@@ -1,6 +1,7 @@
 .vscode/
 .R~
 .m~
+.pyc
 
 # Byte-compiled / optimized / DLL files
 __pycache__/

README.md

@@ -1,10 +1,10 @@
 # pmh-tutorial
 
-This code was downloaded from < https://github.com/compops/pmh-tutorial > and contains the code used to produce the results in the tutorial:
+This code was downloaded from https://github.com/compops/pmh-tutorial and contains the code used to produce the results in the tutorial:
 
 J. Dahlin and T. B. Schön, **Getting started with particle Metropolis-Hastings for inference in nonlinear models**. Pre-print, arXiv:1511:01707, 2017. 
 
-The papers are available as a preprint from < http://arxiv.org/pdf/1511.01707 >.
+The papers are available as a preprint from http://arxiv.org/pdf/1511.01707
 
 Included material
 --------------
@@ -18,4 +18,28 @@ Included material
 
 **matlab-skeleton** Skeleton code files for MATLAB to help step-by-step implementation during courses and seminars. 
 
-
+Copyright information
+--------------
+See *LICENSE* for more information.
+
+``` R
+##############################################################################
+#
+# Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+##############################################################################
+```

matlab/README.md

@@ -1,12 +1,12 @@
 # MATLAB code for PMH tutorial
 
-This MATLAB code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatilty (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at: < http://arxiv.org/pdf/1511.01707 >.
+This MATLAB code implements the Kalman filter (KF), particle filter (PF) and particle Metropolis-Hastings (PMH) algorithm for two different dynamical models: a linear Gaussian state-space (LGSS) model and a stochastic volatilty (SV) model. Note that the Kalman filter can only be employed for the first of these two models. The details of the code is described in the tutorial paper available at http://arxiv.org/pdf/1511.01707
 
-Note that the MATLAB code in this folder covers the basic implementations in the paper. See the R code in r/ for all the implementations and to recreate the results in the tutorial.
+Note that the MATLAB code in this folder covers the basic implementations in the paper. The notation of the variables has been changed sligthly compared with the tutorial paper to improve readability of the code. However, it should be easy to translate between the two. See the R code in r/ for all the implementations and to recreate the results in the tutorial.
 
 Requirements
 --------------
-The code is written and tested for MATLAB 2016b and makes use of the statistics toolbox and the Quandl package. See < https://github.com/quandl/Matlab > for more installation and to download the toolbox. Note that urlread2 is required by the Quandl toolbox and should be installed as detailed in the README file of the Quandl toolbox.
+The code is written and tested for MATLAB 2016b and makes use of the statistics toolbox and the Quandl package. See https://github.com/quandl/Matlab for more installation and to download the toolbox. Note that urlread2 is required by the Quandl toolbox and should be installed as detailed in the README file of the Quandl toolbox.
 
 Main script files
 --------------
@@ -29,4 +29,28 @@ Supporting files
 
 Adapting the code for another model
 --------------
-See the discussion in *README.MD* in the directory *r/*.
+See the discussion in *README.MD* in the directory *r/*.
+
+Copyright information
+--------------
+``` R
+##############################################################################
+#
+# Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
+#
+# This program is free software; you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation; either version 2 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License along
+# with this program; if not, write to the Free Software Foundation, Inc.,
+# 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+#
+##############################################################################
+```

matlab/example1_lgss.m

@@ -1,87 +1,39 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-% Example of state estimation in a LGSS model 
-% using particle filters and Kalman filters
-%
-% Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
-%
-% This program is free software; you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation; either version 2 of the License, or
-% (at your option) any later version.
-%
-% This program is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License along
-% with this program; if not, write to the Free Software Foundation, Inc.,
-% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% State estimation in a LGSS model using particle and Kalman filters
 
 % Set random seed
 rng(0)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Define the model
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Here, we use the following model
-%
-% x[t + 1] = phi * x[t] + sigmav * v[t]
-% y[t] = x[t] + sigmae * e[t]
-%
-% where v[t] ~ N(0, 1) and e[t] ~ N(0, 1)
-
-% Set the parameters of the model
+% x[t + 1] = phi * x[t] + sigmav * v[t],    v[t] ~ N(0, 1)
+% y[t] = x[t] + sigmae * e[t],              e[t] ~ N(0, 1)
 phi = 0.75;
 sigmav = 1.00;
 sigmae = 0.10;
-theta = [phi sigmav sigmae];
-
-% Set the number of time steps to simulate
-T = 250;
-
-% Set the initial state
+parameters = [phi sigmav sigmae];
+noObservations = 250;
 initialState = 0;
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Generate data
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[states, observations] = generateData(parameters, noObservations, initialState);
 
-[x, y] = generateData(theta, T, initialState);
-
-% Plot the measurements and latent states
 subplot(3,1,1); 
-plot(y(2:(T + 1)), 'LineWidth', 1.5, 'Color', [27 158 119] / 256); 
+plot(observations(2:(noObservations + 1)), 'LineWidth', 1.5, 'Color', [27 158 119] / 256); 
 xlabel('time'); 
 ylabel('measurement');
 
 subplot(3,1,2); 
-plot(x(2:(T + 1)), 'LineWidth', 1.5, 'Color', [217 95 2] / 256); 
+plot(states(2:(noObservations + 1)), 'LineWidth', 1.5, 'Color', [217 95 2] / 256); 
 xlabel('time'); 
 ylabel('latent state');
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% State estimation
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Using N = 20 particles and plot the estimate of the latent state
-outPF  = particleFilter(y, theta, 20, initialState);
+% State estimation using particle filter with N = 20 particles
+stateEstPF = particleFilter(observations, parameters, 20, initialState);
 
-% Kalman filter
-outKF  = kalmanFilter(y, theta, initialState, 0.01);
+% State estimation using Kalman filter
+stateEstKF = kalmanFilter(observations, parameters, initialState, 0.01);
 
-% Plot the difference
 subplot(3,1,3); 
-difference = outPF(2:T) - outKF(2:T);
-plot(1:(T - 1), difference, 'LineWidth', 1.5, 'Color', [117 112 179] / 256);
+difference = stateEstPF(2:noObservations) - stateEstKF(2:noObservations);
+plot(1:(noObservations - 1), difference, 'LineWidth', 1.5, 'Color', [117 112 179] / 256);
 xlabel('time'); 
-ylabel('difference in state estimate');
-
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% End of file
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ylabel('difference in state estimate');

matlab/example2_lgss.m

@@ -1,120 +1,62 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
 % Example of particle Metropolis-Hastings in a LGSS model.
-%
-% Copyright (C) 2017 Johan Dahlin < liu (at) johandahlin.com.nospam >
-%
-% This program is free software; you can redistribute it and/or modify
-% it under the terms of the GNU General Public License as published by
-% the Free Software Foundation; either version 2 of the License, or
-% (at your option) any later version.
-%
-% This program is distributed in the hope that it will be useful,
-% but WITHOUT ANY WARRANTY; without even the implied warranty of
-% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-% GNU General Public License for more details.
-%
-% You should have received a copy of the GNU General Public License along
-% with this program; if not, write to the Free Software Foundation, Inc.,
-% 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 % Set random seed
 rng(0)
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Define the model
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-% Here, we use the following model
-%
-% x[t + 1] = phi * x[t] + sigmav * v[t]
-% y[t] = x[t] + sigmae * e[t]
-%
-% where v[t] ~ N(0, 1) and e[t] ~ N(0, 1)
-
-% Set the parameters of the model
+% x[t + 1] = phi * x[t] + sigmav * v[t],    v[t] ~ N(0, 1)
+% y[t] = x[t] + sigmae * e[t],              e[t] ~ N(0, 1)
 phi = 0.75;
 sigmav = 1.00;
 sigmae = 0.10;
-theta = [phi sigmav sigmae];
-
-% Set the number of time steps to simulate
-T = 250;
-
-% Set the initial state
+parameters = [phi sigmav sigmae];
+noObservations = 250;
 initialState = 0;
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Generate data
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
-[x, y] = generateData(theta, T, initialState);
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% Parameter estimation using PMH
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+[states, observations] = generateData(parameters, noObservations, initialState);
 
-% The inital guess of the parameter
+% Setings for PMH
 initialPhi = 0.50;
-
-% No. particles in the particle filter (choose noParticles ~ T)
-noParticles = 100;
-
-% The length of the burn-in and the no. iterations of PMH 
-% (noBurnInIterations < noIterations)
+noParticles = 100;          % Use noParticles ~ noObservations
 noBurnInIterations = 1000;
 noIterations = 5000;
-
-% The standard deviation in the random walk proposal
 stepSize = 0.10;
 
 % Run the PMH algorithm
-res = particleMetropolisHastings(y, initialPhi, [sigmav sigmae],...
-                                 noParticles, initialState,...
-                                 noIterations, stepSize);
+phiTrace = particleMetropolisHastings(observations, initialPhi, [sigmav sigmae], noParticles, initialState, noIterations, stepSize);
 
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % Plot the results
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-
 noBins = floor(sqrt(noIterations - noBurnInIterations));
 grid = noBurnInIterations:noIterations;
-resPhi = res(noBurnInIterations:noIterations);
+phiTrace = phiTrace(noBurnInIterations:noIterations);
 
-% Plot the parameter posterior estimate
-% Solid black line indicate posterior mean
+% Plot the parameter posterior estimate (solid black line = posterior mean)
 subplot(3, 1, 1);
-hist(resPhi, noBins);
+hist(phiTrace, noBins);
 xlabel('phi'); 
 ylabel('posterior density estimate');
 
 h = findobj(gca, 'Type', 'patch');
 set(h, 'FaceColor', [117 112 179] / 256, 'EdgeColor', 'w');
 
 hold on;
-plot([1 1] * mean(resPhi), [0 200], 'LineWidth', 3);
+plot([1 1] * mean(phiTrace), [0 200], 'LineWidth', 3);
 hold off;
 
-% Plot the trace of the Markov chain after burn-in
-% Solid black line indicate posterior mean
+% Plot the trace of the Markov chain after burn-in  (solid black line = posterior mean)
 subplot(3, 1, 2);
-plot(grid, resPhi, 'Color', [117 112 179] / 256, 'LineWidth', 1);
+plot(grid, phiTrace, 'Color', [117 112 179] / 256, 'LineWidth', 1);
 xlabel('iteration'); 
 ylabel('phi');
 
 hold on;
-plot([grid(1) grid(end)], [1 1] * mean(resPhi), 'k', 'LineWidth', 3);
+plot([grid(1) grid(end)], [1 1] * mean(phiTrace), 'k', 'LineWidth', 3);
 hold off;
 
 % Plot ACF of the Markov chain after burn-in
 subplot(3, 1, 3);
-[acf, lags] = xcorr(resPhi - mean(resPhi), 100, 'coeff');
+[acf, lags] = xcorr(phiTrace - mean(phiTrace), 100, 'coeff');
 stem(lags(101:200), acf(101:200), 'Color', [117 112 179] / 256, 'LineWidth', 2);
 xlabel('lag'); 
-ylabel('ACF of phi');
-
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-% End of file
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ylabel('ACF of phi');