octave_packages/optim-1.2.0/LinearRegression.m

   1 ## Copyright (C) 2007 Andreas Stahel <Andreas.Stahel@bfh.ch>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## general linear regression
  17 ##
  18 ## [p,y_var,r,p_var]=LinearRegression(F,y)
  19 ## [p,y_var,r,p_var]=LinearRegression(F,y,weight)
  20 ##
  21 ## determine the parameters p_j  (j=1,2,...,m) such that the function
  22 ## f(x) = sum_(i=1,...,m) p_j*f_j(x) fits as good as possible to the
  23 ## given values y_i = f(x_i)
  24 ##
  25 ## parameters
  26 ## F  n*m matrix with the values of the basis functions at the support points
  27 ##    in column j give the values of f_j at the points x_i  (i=1,2,...,n)
  28 ## y  n column vector of given values
  29 ## weight  n column vector of given weights
  30 ##
  31 ## return values
  32 ## p     m vector with the estimated values of the parameters
  33 ## y_var estimated variance of the error
  34 ## r     weighted norm of residual
  35 ## p_var estimated variance of the parameters p_j
  36
  37 function [p,y_var,r,p_var]=LinearRegression(F,y,weight)
  38
  39 if (nargin < 2 || nargin >= 4)
  40  usage('wrong number of arguments in [p,y_var,r,p_var]=LinearRegression(F,y)');
  41 end
  42
  43 [rF, cF] = size(F);  [ry, cy] =size(y);
  44 if (rF ~= ry || cy > 1)
  45   error ('LinearRegression: incorrect matrix dimensions');
  46 end
  47
  48 if (nargin==2)  % set uniform weights if not provided
  49   weight=ones(size(y));
  50 end
  51
  52 %% Fw=diag(weight)*F;
  53 wF=F;
  54 for j=1:cF
  55   wF(:,j)=weight.*F(:,j);
  56 end
  57
  58 [Q,R]=qr(wF,0);                  % estimate the values of the parameters
  59 p=R\(Q'*(weight.*y));
  60
  61
  62 residual=F*p-y;                  % compute the residual vector
  63 r=norm(weight.*residual);   % and its weighted norm
  64                                  % variance of the weighted y-errors
  65 y_var= sum((residual.^2).*(weight.^4))/(rF-cF);
  66
  67 if nargout>3    % compute variance of parameters only if needed
  68 %%  M=inv(R)*Q'*diag(weight);
  69   M=inv(R)*Q';
  70   for j=1:cF
  71     M(j,:)=M(j,:).*(weight');
  72   end
  73   M=M.*M;                        % square each entry in the matrix M
  74   p_var=M*(y_var./(weight.^4));  % variance of the parameters
  75 end
  76