updates for the arXiv-version

Former-commit-id: 64843fd [formerly 64843fd [formerly 2debf13]]
Former-commit-id: 4b8afdcec2a00931e4714c18e8d52328225f6d76
Former-commit-id: 0a15fe9

Author: Johan Dahlin

README.md

@@ -2,9 +2,9 @@
 
 This code was downloaded from < https://github.com/compops/pmh-tutorial > or from < http://users.isy.liu.se/en/rt/johda87/ > 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:1510:xxxxv1, 2015. 
+* J. Dahlin and T. B. Schön, *Getting started with particle Metropolis-Hastings for inference in nonlinear models*. Pre-print, arXiv:1511:xxxxv1, 2015. 
 
-The papers are available as a preprint from < http://arxiv.org/pdf/1510.xxxxv1 > and < http://users.isy.liu.se/en/rt/johda87/ >.
+The papers are available as a preprint from < http://arxiv.org/pdf/1511.xxxxv1 > and < http://users.isy.liu.se/en/rt/johda87/ >.
 
 Requirements
 --------------
@@ -22,14 +22,14 @@ Included files
 
 **example4-sv.[R,py,m]** Implements the numerical illustration in Section 4.2, which makes use of the same setup as in Section 4.1 but with a tuned proposal distribution. The output is the estimated ACF and IACT as presented in Figure 7. The code takes some time (hours to execute).
 
-**example*.[RData,mat]** A saved copy of the workspace after running to code.
+**example*.[RData,mat,spdata]** A saved copy of the workspace after running thr corresponding code. Can be used to directly recreate the plots in the tutorial and to conduct additional analysis.
 
 Supporting files
 --------------
-**stateEstimationHelper.[py,m]**
+**stateEstimationHelper.[py,R]**
 Implementes the data generation for the LGSS model (generateData), the faPF for the LGSS model (sm), the Kalman filter for the LGSS model (kf) and the bPF for the SV model (sm_sv). In Matlab, these functions are defined in four seperate m-files with the corresponding file names.
 
-**parameterEstimationHelper.[py,m]**
+**parameterEstimationHelper.[py,R]**
 Implementes the PMH algorithm for the LGSS model (pmh) and the SV model (pmh_sv). In Matlab, these functions are defined in two seperate m-files with the corresponding file names.
 
 **omxs30data.csv**

matlab/example1_lgss.mat

@@ -65,7 +65,7 @@
 % Solid black line indicate posterior mean
 subplot(4,2,3);
 hist(thhat(nBurnIn:nRuns,1), floor(sqrt(nRuns-nBurnIn)) );
-xlabel('phi'); ylabel('posterior density estimate');
+xlabel('mu'); ylabel('posterior density estimate');
 
 h = findobj(gca,'Type','patch');
 set(h,'FaceColor',[117 112 179]/256,'EdgeColor','w');
@@ -74,24 +74,24 @@
 subplot(4,2,4);
 plot(nBurnIn:nRuns,thhat(nBurnIn:nRuns,1), 'Color', [117 112 179]/256);
 hold on; plot([nBurnIn nRuns],[1 1] * mean(thhat(nBurnIn:nRuns,1)), 'k'); hold off;
-xlabel('iteration'); ylabel('trace of phi');
+xlabel('iteration'); ylabel('trace of mu');
 
 subplot(4,2,5);
 hist(thhat(nBurnIn:nRuns,2), floor(sqrt(nRuns-nBurnIn)) );
-xlabel('sigmav'); ylabel('posterior density estimate');
+xlabel('phi'); ylabel('posterior density estimate');
 
 h = findobj(gca,'Type','patch');
 set(h,'FaceColor',[231 41 138]/256,'EdgeColor','w');
 hold on; plot([1 1] * mean(thhat(nBurnIn:nRuns,2)), [0 500], 'k'); hold off;
 
 subplot(4,2,6);
 plot(nBurnIn:nRuns,thhat(nBurnIn:nRuns,2), 'Color', [231 41 138]/256);
-xlabel('iteration'); ylabel('trace of sigmav');
+xlabel('iteration'); ylabel('trace of phi');
 hold on; plot([nBurnIn nRuns],[1 1] * mean(thhat(nBurnIn:nRuns,2)), 'k'); hold off;
 
 subplot(4,2,7);
 hist(thhat(nBurnIn:nRuns,3), floor(sqrt(nRuns-nBurnIn)) );
-xlabel('beta'); ylabel('posterior density estimate');
+xlabel('sigmav'); ylabel('posterior density estimate');
 
 h = findobj(gca,'Type','patch');
 set(h,'FaceColor',[102 166 30]/256,'EdgeColor','w');
@@ -100,7 +100,7 @@
 subplot(4,2,8);
 plot(nBurnIn:nRuns,thhat(nBurnIn:nRuns,3), 'Color', [102 166 30]/256);
 hold on; plot([nBurnIn nRuns],[1 1] * mean(thhat(nBurnIn:nRuns,3)), 'k'); hold off;
-xlabel('iteration'); ylabel('trace of beta');
+xlabel('iteration'); ylabel('trace of sigmav');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % End of file

matlab/example2_lgss.mat

@@ -0,0 +1 @@
+d5935fe8d440d8277110a00f97551bf398434497

matlab/example3_sv.m

@@ -69,7 +69,7 @@
 % Solid black line indicate posterior mean
 subplot(4,2,3);
 hist(thhat(nBurnIn:nRuns,1), floor(sqrt(nRuns-nBurnIn)) );
-xlabel('phi'); ylabel('posterior density estimate');
+xlabel('mu'); ylabel('posterior density estimate');
 
 h = findobj(gca,'Type','patch');
 set(h,'FaceColor',[117 112 179]/256,'EdgeColor','w');
@@ -78,24 +78,24 @@
 subplot(4,2,4);
 plot(nBurnIn:nRuns,thhat(nBurnIn:nRuns,1), 'Color', [117 112 179]/256);
 hold on; plot([nBurnIn nRuns],[1 1] * mean(thhat(nBurnIn:nRuns,1)), 'k'); hold off;
-xlabel('iteration'); ylabel('trace of phi');
+xlabel('iteration'); ylabel('trace of mu');
 
 subplot(4,2,5);
 hist(thhat(nBurnIn:nRuns,2), floor(sqrt(nRuns-nBurnIn)) );
-xlabel('sigmav'); ylabel('posterior density estimate');
+xlabel('phi'); ylabel('posterior density estimate');
 
 h = findobj(gca,'Type','patch');
 set(h,'FaceColor',[231 41 138]/256,'EdgeColor','w');
 hold on; plot([1 1] * mean(thhat(nBurnIn:nRuns,2)), [0 500], 'k'); hold off;
 
 subplot(4,2,6);
 plot(nBurnIn:nRuns,thhat(nBurnIn:nRuns,2), 'Color', [231 41 138]/256);
-xlabel('iteration'); ylabel('trace of sigmav');
+xlabel('iteration'); ylabel('trace of phi');
 hold on; plot([nBurnIn nRuns],[1 1] * mean(thhat(nBurnIn:nRuns,2)), 'k'); hold off;
 
 subplot(4,2,7);
 hist(thhat(nBurnIn:nRuns,3), floor(sqrt(nRuns-nBurnIn)) );
-xlabel('beta'); ylabel('posterior density estimate');
+xlabel('sigmav'); ylabel('posterior density estimate');
 
 h = findobj(gca,'Type','patch');
 set(h,'FaceColor',[102 166 30]/256,'EdgeColor','w');
@@ -104,7 +104,7 @@
 subplot(4,2,8);
 plot(nBurnIn:nRuns,thhat(nBurnIn:nRuns,3), 'Color', [102 166 30]/256);
 hold on; plot([nBurnIn nRuns],[1 1] * mean(thhat(nBurnIn:nRuns,3)), 'k'); hold off;
-xlabel('iteration'); ylabel('trace of beta');
+xlabel('iteration'); ylabel('trace of sigmav');
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 % End of file

matlab/example3_sv.mat.REMOVED.git-id

@@ -0,0 +1 @@
+eecf624ca906846eea8c070b8b05254e79899f0b

matlab/example4_sv.m

@@ -29,13 +29,14 @@
 function[xhatf,llp] = sm_sv(y,mu,phi,sigmav,N,T)
 
   % Initalise variables
-  xhatf = zeros( T+1, 1);
+  a     = zeros( N, T+1);
   p     = zeros( N, T+1);
+  v     = ones(  N, T+1);
   w     = ones(  N, T+1) / N;
   llp   = 0;
 
+  a(:,1)   = 1:N;
   p(:,1)   = mu + sigmav/sqrt(1-phi^2) * normrnd( 0, 1, N, 1 ) ;
-  xhatf(1) = mean(p(:,1));
 
   %===========================================================
   % Run main loop
@@ -48,6 +49,12 @@
     idx = randsample( N, N, true, w(:,tt-1) ); 
     idx = idx( randperm( N ) );
 
+    % Resample the ancestory linage
+    a(:,1:tt-1) = a(idx,1:tt-1);
+    
+    % Add the most recent ancestors
+    a(:,tt)     = idx;    
+    
     %=========================================================
     % Propagate
     %=========================================================
@@ -56,23 +63,24 @@
     %=========================================================
     % Compute weights
     %=========================================================
-    w(:,tt) = dnorm( y(tt-1), 0, exp(p(:,tt)/2) );
+    v(:,tt) = dnorm( y(tt-1), 0, exp(p(:,tt)/2) );
 
     % Rescale log-weights and recover weight
-    wmax    = max( w(:,tt) );
-    w(:,tt) = exp( w(:,tt) - wmax );
+    vmax    = max( v(:,tt) );
+    v(:,tt) = exp( v(:,tt) - vmax );
 
     % Estimate the log-likelihood
-    llp    = llp + wmax + log(sum( w(:,tt) )) - log(N);
+    llp    = llp + vmax + log(sum( v(:,tt) )) - log(N);
 
     % Normalize the weights
-    w(:,tt) = w(:,tt) / sum( w(:,tt) );
-    
-    % Estimate the state
-    xhatf(tt) = sum( w(:,tt) .* p(:,tt) );
+    w(:,tt) = v(:,tt) / sum( v(:,tt) );
 
   end
 
+  % Sample the state estimate using the weights at tt=T
+  nIdx  = randsample( N, 1, true, w(:,T) );
+  indices = sub2ind(size(p), a(nIdx,:), 1:(T+1));
+  xhatf = p( indices );
 end
 
 % Helper for computing the logarithm of the Gaussian density