//----------------------------------------------------------------------------//
// linfit: linear mutivariate regression for Scilab http://www.scilab.org/    //
// Author: Javier I. Carrero, jicarrerom@unal.edu.co                          //
// Last modification date: 2012-10-5                                          //
// Warning this function is provided as-is, the author does not provide any   //
// guarantee about its results.                                               //
//----------------------------------------------------------------------------//
// Mathematical background:                                                   //
// Based on the minimization of residuals, see polyfit documentation          //
// and Chapra and Canale's "Numerical Methods for Engineers", 5th ed., ch. 17 //
// (McGraw-Hill, 2005) ISBN-13: 978-0073101569                                //
//----------------------------------------------------------------------------//
// Default behavior:                                                          //
// For y that depends on n independent variables x1, x2, ..., xn this         //
// function calculates the values of the coefficients ai that best fit        //
// y = a0 + (a1*x1) + (a2*x2) + ... + (an*xn )                                //
// in the sense of minimizing the sum of square residuals calculated from a   //
// data set of m elements (y; x1, x2, ...xn)                                  //
// Alternate behavior:                                                        //
// When invoked with the noA0 option linfit forces a0 = 0 producing           //
// y = (a1*x1) + (a2*x2) + ... + (an*xn )                                     //
//----------------------------------------------------------------------------//
// Arguments (function input)                                                 //
// Given m (y; x1, x2, ...xn) sets of data for n independent x variables:     //
// * xVar is a (m, n)  matrix. Each column in xVar contains the m values      //
//   their respective variable. For example xVar(:,1) = [x11, x21, ... xm1]'  //
// * yVar is a (m, 1) matrix, i.e. a column vector. Each element represents   //
//   a y value tied to the respective column in xVar. Given yVar = [y1 y2 ... //
//   ym]'  y1 corresponds to xVar(1,:), y2 to xVar(2,:), and so on.           //
// * noA0 (optional): if present produces an output with a0=0. To invoke this //
//   option add this argument: noA0 = "Y"                                     //
// * doGraph (optional): if present produces a graphic in which the           //
//   calculated values of y are plotted against the perfect fit, i.e.         //
//   against the yCalc=y line. To invoke this option add this argument:       //
//   doGraph = "yVar"                                                         // 
//----------------------------------------------------------------------------//
// Results (function output):                                                 //
// * listCoefs in the default form is a (n+1, 1) matrix containing the values //
//   of ai, i.e. listCoefs = [a0 a1 a2 ... an]'. If linfit is invoked with    //
//   the noA0 option listCoefs will contain n values, i.e.                    //
//   listCoefs = [a1 a2 ... an]'                                              //
// * yCalc is a column vector with m elements corresponding to the y values   //
//   calculated with the m sets of x1, x2, ..., xn                            //
// * statParams is a vector with statistical parameters defined as            //
//   statParams = [St Sr stdv r2] where                                       //
//   - St: sum of the square differences between yi and yavg, being yavg      //
//     the average of y                                                       //
//   - Sr: sum of the m residuals, being each residual (yi-yi_calc)^2         //
//     being yCalc the value obtained from the set x1i, x2i, ..., xni         //
//   - stdv: standard deviation, defined as ( St / m-1 ) ^ 0.5                //
//   - r2: correlation coefficient, defined as r2 = ( St - Sr ) / St          //
//   - Syx: standard error, defined as Syx = (Sr/(m-(n+1)))^1/2               //
//----------------------------------------------------------------------------//

function [listCoefs, yCalc, statParams] = linfit(xVar, yVar, noA0, doGraph)

//----------------------------------------------------------------------------//
// Argument Checking                                                          //
//----------------------------------------------------------------------------//
[lhs, rhs] = argn()

if ( ( rhs < 2 ) | ( rhs > 4 ) ) then
  error(" Wrong number of input arguments. ")
end

if ( lhs > 3 ) then
  error(" Wrong number of output arguments. ")
end

if ( typeof(xVar) <> "constant") then
  error("xVar should be a vector of real numbers")
end

if ( typeof(yVar) <> "constant") then
  error("yVar should be a vector of real numbers")
end

//----------------------------------------------------------------------------//
// Note:                                                                      // 
// numPoints = m  (number of data points)                                     // 
// numVar = n     (number of independent variables                            //
//----------------------------------------------------------------------------//
[numPoints, numVar]  = size(xVar)
[numRowsY, numColsY] = size(yVar)

if ( numRowsY <> numPoints ) then
  error(" Number of rows of yVar must be equal to the number of columns of xVar")
end,

if ( numColsY <> 1 ) then
  error(" Inconsistent argument: yVar must be a column vector.")
end,

if ( numPoints < ( numVar + 1 ) ) then
   printf(" There are less data points than variables+1, check")
end  

// check for the optional parameters
if exists("noA0") then
  calculateA0 = %f
else
  calculateA0 = %t
end,

if exists("doGraph") then
  doPlot = %t
else
  doPlot = %f
end,

//----------------------------------------------------------------------------//
// Building of Z matrix                                                       //
// If calculateA0== %f the alternate version (a0=0) is produced               //
//----------------------------------------------------------------------------//
if ( calculateA0 ) then
  Z = zeros(numPoints, numVar+1)
  Z(:,1) = 1
  Z(:,2:numVar+1) = xVar
else
  Z = zeros(numPoints, numVar)
  Z = xVar
end,

//----------------------------------------------------------------------------//
// Main output                                                                //
//----------------------------------------------------------------------------//
// The coefficients (ai) are grouped in the variable listCoefs, which         //
// correspond to the result                                                   //
// listCoefs = ( Z' * Z ) \ ( Z' * yVar )                                     //
// but given that in some cases Z, yVar produce ill conditioned               //
// matrices it was decided to use the Q-R method of solution                  //
// It takes a longer time (not very much) but is reliable                     //
//----------------------------------------------------------------------------//

izq = Z' * Z
der = Z' * yVar

[Q, R] = qr(izq)

Qtb = Q' * der

if ( calculateA0 ) then
    n = numVar + 1
else
    n = numVar
end,

listCoefs    = zeros(n,1)
listCoefs(n) = Qtb(n) / R(n,n)

for j=n-1:-1:1
    listCoefs(j) = ( Qtb(j) - ( R(j,j+1:n) * listCoefs(j+1:n) ) ) / R(j,j)
end,

yCalc = zeros(numPoints,1)   

if ( calculateA0 ) then
    for i=1:numPoints
        yCalc(i) = sum( listCoefs' .* [1 xVar(i,:)] )
    end,  
else
    for i=1:numPoints
        yCalc(i) = sum( listCoefs' .* xVar(i,:) )
    end,    
end,

//----------------------------------------------------------------------------//
// Statistics                                                                 //
//----------------------------------------------------------------------------//
avg_y = sum( yVar ) / numPoints
St    = sum( ( yVar - avg_y ) .^ 2 )
Sr    = sum( ( yCalc - yVar ) .^ 2 )
stdv  = sqrt( St / ( numPoints - 1 ) )

if ( St >= %eps ) then
      r2 = ( St - Sr ) / St
else
    r2 = %inf
end

if ( numPoints > ( numVar + 1 ) ) then
    Syx = sqrt( Sr / ( numPoints - (numVar + 1 ) ) )
else
    Syx = %inf
end,

statParams = [St Sr stdv r2 Syx]';

//----------------------------------------------------------------------------//
// Graphic checking (optional, if  doGraph is present)                        //
//----------------------------------------------------------------------------//
if (doPlot) then
    numTestPoints = 100
    ytest = linspace(min(yVar), max(yVar), numTestPoints)

    clf()
    plot(yVar, yCalc, 'o', ytest, ytest, '-')
    ejes = get("current_axes")
    ejes.x_label.text = 'y'
    ejes.y_label.text = 'y_calc'
end,
endfunction
