Added function that does the equivalent of orthogonal projection for …

…kernels, for some ideas on nonlinear sysid I want to try.

glm/gen_glm.m

@@ -1,19 +1,20 @@
 function [x y] = gen_glm(A,B,C,f,b,x0,u,Q)
 
 n = size(x0,1);
-x = zeros(n,size(u,2));
+T = size(u,2);
+x = zeros(n,T);
 
 if isempty(Q)
     Q = zeros(n);
 end
 
 x(:,1) = x0;
 E = chol(Q);
-for i = 1:size(u,2)
+for i = 1:T
     e = E*randn(n,1);
     if i < size(u,2)
         x(:,i+1) = A*x(:,i) + B*u(:,i) + e(1:n);
     end
 end
 
-y = poissrnd( f( C*x + b ) );
+y = poissrnd( f( C*x + b*ones(1,T) ) );

glm/glm_ekf.m

@@ -26,7 +26,7 @@
 for i = 1:T
     R = f(C*zt + b);
     Rinv = diag(1./R);
-    H = grad(C*xt + b)*ones(1,k) .* C;
+    H = grad(C*zt + b)*ones(1,k) .* C;
     % Precision matrix of y(:,i) given y(:,1:i-1).
     if size(y,1) > 50
         % Same as (H*Pt*H' + R)^-1 by Woodbury lemma.

kernel_ssid/orth_kernel.m

@@ -0,0 +1,19 @@
+function ku = orth_kernel(k,U)
+% given a kernel k(x,x') and a matrix U, construct a kernel ku that is the
+% inner product in the orthogonal projection of U in feature space.  That
+% is, ku(x,u) = 0 for all x and u such that u is a column of U.
+
+Kuu = zeros(size(U,2));
+for i = 1:size(U,2)
+    for j = 1:size(U,2)
+        Kuu(i,j) = k(U(:,i),U(:,j));
+    end
+end
+invKuu = Kuu^-1;
+ku = @(x,y) k(x,y) - kux(x,U,k)'*invKuu*kux(y,U,k);
+
+function z = kux(x,U,k)
+z = zeros(size(U,2));
+for i = 1:size(U,2)
+    z(i) = k(x,U(:,i));
+end