From 6efa57a692372e0bf91ea2b389ca38f3750aa37d Mon Sep 17 00:00:00 2001
From: Michel Juillard <michel.juillard@ens.fr>
Date: Mon, 22 Feb 2010 09:29:47 +0100
Subject: [PATCH] corrected bug for decompostion (nargout > 7) in
 missing_DiffuseKalmanSmoother3.m (cherry picked from commit
 f4b7840739743f15fa763de1c527277492dbfc73)

---
 matlab/missing_DiffuseKalmanSmoother3.m | 16 +++++++++++-----
 1 file changed, 11 insertions(+), 5 deletions(-)

diff --git a/matlab/missing_DiffuseKalmanSmoother3.m b/matlab/missing_DiffuseKalmanSmoother3.m
index 8f2a04151e..3ec17d6a50 100644
--- a/matlab/missing_DiffuseKalmanSmoother3.m
+++ b/matlab/missing_DiffuseKalmanSmoother3.m
@@ -292,19 +292,25 @@ if nargout > 7
     decomp = zeros(nk,mm,rr,smpl+nk);
     ZRQinv = inv(Z*QQ*Z');
     for t = max(d,1):smpl
-        ri_d = Z'*iF(:,:,t)*v(:,t);
+        ri_d = zeros(mm,1);
+        di = flipud(data_index{t})';
+        for i = di
+            if Fi(i,t) > crit
+                ri_d = Z(i,:)'/Fi(i,t)*v(i,t)+ri_d-Ki(:,i,t)'*ri_d/Fi(i,t)*Z(i,:)';
+            end
+        end
         
         % calculate eta_tm1t
         eta_tm1t = QRt*ri_d;
         % calculate decomposition
         Ttok = eye(mm,mm); 
+        AAA = P1(:,:,t)*Z'*ZRQinv*Z*R;
         for h = 1:nk
+            BBB = Ttok*AAA;
             for j=1:rr
-                eta=zeros(rr,1);
-                eta(j) = eta_tm1t(j);
-                decomp(h,:,j,t+h) = T^(h-1)*P(:,:,t)*Z'*ZRQinv*Z*R*eta;
+                decomp(h,:,j,t+h) = eta_tm1t(j)*BBB(:,j);
             end
+            Ttok = T*Ttok;
         end
     end
 end
-
-- 
GitLab