making univariate Kalman filter code simpler and more efficient

matlab/kalman/likelihood/univariate_kalman_filter.m

@@ -132,27 +132,25 @@ while notsteady && t<=last
     for i=1:rows(z)
         if Zflag
             prediction_error = Y(d_index(i),t) - z(i,:)*a;
-            Fi = z(i,:)*P*z(i,:)' + H(d_index(i));
+            PZ = P*z(i,:)';
+            Fi = z(i,:)*PZ + H(d_index(i));
         else
             prediction_error = Y(d_index(i),t) - a(z(i));
-            Fi = P(z(i),z(i)) + H(d_index(i));
+            PZ = P(:,z(i));
+            Fi = PZ(z(i)) + H(d_index(i));
         end
         if Fi>kalman_tol
-            if Zflag
-                Ki =  (P*z(i,:)')/Fi;
-            else
-                Ki = P(:,z(i))/Fi;
-            end
+            Ki =  PZ/Fi;
             if t>no_more_missing_observations
                 K(:,i) = Ki;
             end
             a = a + Ki*prediction_error;
-            P = P - (Fi*Ki)*Ki';
+            P = P - PZ*Ki';
             lik(s) = lik(s) + log(Fi) + prediction_error*prediction_error/Fi + l2pi;
         end
     end
     a = T*a;
-    P = T*P*transpose(T) + QQ;
+    P = T*P*T' + QQ;
     if t>=no_more_missing_observations
         notsteady = max(abs(K(:)-oldK))>riccati_tol;
         oldK = K(:);

matlab/kalman/likelihood/univariate_kalman_filter_ss.m

@@ -90,19 +90,17 @@ while t<=last
     for i=1:pp
         if Zflag
             prediction_error = Y(i,t) - Z(i,:)*a;
-            Fi = Z(i,:)*PP*Z(i,:)' + H(i);
+            PPZ = PP*Z(i,:)';
+            Fi = Z(i,:)*PPZ + H(i);
         else
             prediction_error = Y(i,t) - a(Z(i));
-            Fi = PP(Z(i),Z(i)) + H(i);
+            PPZ = PP(:,Z(i));
+            Fi = PPZ(Z(i)) + H(i);
         end
         if Fi>kalman_tol
-            if Zflag
-                Ki = (PP*Z(i,:))'/Fi;
-            else
-                Ki = PP(:,Z(i))/Fi;
-            end
+            Ki = PPZ/Fi;
             a  = a + Ki*prediction_error;
-            PP = PP - (Fi*Ki)*transpose(Ki);
+            PP = PP - PPZ*Ki';
             likk(s) = likk(s) + log(Fi) + prediction_error*prediction_error/Fi + l2pi;
         end
     end