projectionmatrix

 function f = projectionmatrix(X)

 <X> is samples x parameters

 what we want to do is to perform a regression using <X>
 and subtract out the fit.  this is accomplished by
 y-X*inv(X'*X)*X'*y = (I-X*inv(X'*X)*X')*y = f*y
 where y is the data (samples x cases).

 what this function does is to return <f> which has
 dimensions samples x samples.  to accomplish this,
 we rely heavily on olsmatrix.m.

 if <X> has no parameters, the output of this function is 1.

 history:
 - 2013/08/18 - in the cases of empty <X>, we now return 1 instead of [].

 example:
 x = sort(randn(100,1));
 x2 = projectionmatrix(constructpolynomialmatrix(100,0:1))*x;
 figure; hold on; plot(x,'r-'); plot(x2,'g-');

0001 function f = projectionmatrix(X)
0002 
0003 % function f = projectionmatrix(X)
0004 %
0005 % <X> is samples x parameters
0006 %
0007 % what we want to do is to perform a regression using <X>
0008 % and subtract out the fit.  this is accomplished by
0009 % y-X*inv(X'*X)*X'*y = (I-X*inv(X'*X)*X')*y = f*y
0010 % where y is the data (samples x cases).
0011 %
0012 % what this function does is to return <f> which has
0013 % dimensions samples x samples.  to accomplish this,
0014 % we rely heavily on olsmatrix.m.
0015 %
0016 % if <X> has no parameters, the output of this function is 1.
0017 %
0018 % history:
0019 % - 2013/08/18 - in the cases of empty <X>, we now return 1 instead of [].
0020 %
0021 % example:
0022 % x = sort(randn(100,1));
0023 % x2 = projectionmatrix(constructpolynomialmatrix(100,0:1))*x;
0024 % figure; hold on; plot(x,'r-'); plot(x2,'g-');
0025 
0026 % handle corner case
0027 if isempty(X)
0028   f = 1;
0029   return;
0030 end
0031 
0032 % do it
0033 f = eye(size(X,1)) - X*olsmatrix(X);