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);