share/lsquares/plsquares.mac

   1 /*
   2 plsquares.mac v1.1 for Maxima (tested with Maxima 5.9.0).
   3
   4 Multivariable polynomial adjustment of a data table by the "least squares"
   5 method.
   6
   7 Usage:
   8   plsquares(Mat, VarList, depvars, maxexpon, maxdegree);
   9     Mat       - a matrix containing the data.
  10     VarList   - list of variable names (one for each Mat column).
  11                 Use "-" instead of varnames to ignore Mat columns.
  12     depvars   - the name of a dependent variable or a list with one or more
  13                 names of dependent variables. The names must be in VarList.
  14     maxexpon  - optional maximum exponent for each independent variable.
  15                 Default: 1.
  16     maxdegree - optional maximum polynomial degree (the sum of exponents of
  17                 each term will be equal or smaller than maxdegree).
  18                 If maxdgree = 0 then no limit is applied.
  19                 Default: maxexpon.
  20
  21 Examples:
  22   The file "plsquares.dem" shows some usage examples.
  23
  24 Results:
  25   - If depvars is the name of a dependent variable (not in a list),
  26   plsquares returns the adjusted polynomial.
  27     If depvars is a list of one or more dependent variables, plsquares
  28   returns a list with the adjusted polynomial(s).
  29   - The Determination Coefficients are displayed in order to inform about
  30   the adjustment goodness (from 0:no correlation to 1:exact correlation).
  31   These values are also stored in the global variable DETCOEF (a list if
  32   depvars is a list).
  33
  34 Dependences:
  35   makeOrders.mac
  36
  37 History:
  38   2003-11 Salvador Bosch Pérez - version 1.1. Multiple dependent variables
  39     (to return a list of polynomials). maxexpon and maxdegree are now
  40     optional. Code more readable.
  41   2003-10 Salvador Bosch Pérez - version 1.0 (not released)
  42
  43 Possible future improvements:
  44   - Option to read the data from a file instead of from a matrix.
  45   - Option to include a column with rows weights.
  46
  47 --
  48
  49 Copyright (C) 2003  Salvador Bosch Pérez
  50
  51 This library is free software; you can redistribute it and/or
  52 modify it under the terms of the GNU Lesser General Public
  53 License as published by the Free Software Foundation; either
  54 version 2.1 of the License, or (at your option) any later version.
  55
  56 This library is distributed in the hope that it will be useful,
  57 but WITHOUT ANY WARRANTY; without even the implied warranty of
  58 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  59 Lesser General Public License for more details.
  60
  61 You should have received a copy of the GNU Lesser General Public
  62 License along with this library; if not, write to the Free Software
  63 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA
  64 */
  65
  66 load("makeOrders")$   /* to obtain all the combinations of polynomial powers */
  67
  68 plsquares([ArgList]):=
  69   block([nargs, Mat, VarList, depvars, maxexpon, maxdegree,
  70          ndat, nvar, DepVarList, depvarncol, DepCols, PowersList, ncoef,
  71          TransformedData, LinearSystem, DepSums, PolyCoef, PolyList,
  72          idat, ivar, jvar, irow, icol, icoef, degreelimit],
  73
  74     /* Elaboration and depuration of the function arguments */
  75     narg:length(ArgList),
  76     if narg < 3 or narg > 5 then (
  77       print("plsquares: bad number of function arguments (it's required 3, 4 or 5 arguments)."),
  78       return(false)
  79     ),
  80     Mat:ArgList[1],
  81     VarList:ArgList[2],
  82     depvars:ArgList[3],
  83     if narg > 3 then maxexpon:ArgList[4]
  84                 else maxexpon:1,
  85     if narg = 5 then maxdegree:ArgList[5]
  86                 else maxdegree:maxexpon,
  87     ndat:length(Mat),
  88     nvar:length(Mat[1]),
  89     if atom(depvars) then DepVarList:[depvars]
  90                      else DepVarList:depvars,
  91     ndepvar:length(DepVarList),
  92     if length(VarList) # nvar then (
  93       print("plsquares: incorrect number of variable names (", nvar,
  94             "matrix columns but", length(VarList), "variable names)."),
  95       return(false)
  96     ),
  97     for ivar:ndepvar thru 1 step -1 do
  98       if member(DepVarList[ivar], VarList) = false then (
  99         print("plsquares: dependent variable", DepVarList[ivar], "isn't in",
 100               VarList, "."),
 101         DepVarList:delete(DepVarList[ivar], DepVarList),
 102         ndepvar = ndepvar - 1
 103       ),
 104     if ndepvar < 1 then (
 105       print("plsquares: no dependent variables."),
 106       return(false)
 107     ),
 108     if maxexpon < 1 then (
 109       print("plsquares: the maximum variable exponent must be greater than 0."),
 110       return(false)
 111     ),
 112     if maxdegree # 0 and maxdegree < maxexpon then (
 113       print("plsquares: the maximum degree of the polynomial must not be smaller than",
 114             maxexpon),
 115       return(false)
 116     ),
 117     for ivar:nvar thru 1 step -1 do
 118       if VarList[ivar] = "-" then (
 119         Mat:submatrix(Mat, ivar),
 120         nvar:nvar - 1
 121       ),
 122     VarList:delete("-", VarList),
 123     for ivar:1 thru ndepvar do (
 124       depvarncol:ev(for jvar:1 thru nvar do
 125                       if VarList[jvar] = DepVarList[ivar] then return(jvar)),
 126       DepCols[ivar]:col(Mat, depvarncol),
 127       VarList:delete(DepVarList[ivar], VarList),
 128       Mat:submatrix(Mat, depvarncol),
 129       nvar:nvar - 1
 130     ),
 131     PowersList:makeOrders(VarList, makelist(maxexpon, i, 1, nvar)),
 132     if maxdegree > 0 then (
 133       degreelimit(l):=lsum(i,i,l)<=maxdegree,
 134       PowersList:sublist(PowersList, degreelimit)
 135     ),
 136     ncoef:length(PowersList),
 137     if ndat < ncoef then (
 138       print("plsquares: insufficient number of data rows (at least", ncoef,
 139             "are required)."),
 140       return(false)
 141     ),
 142     apply(kill, VarList),
 143
 144     /* Preparation of the linear system */
 145     LinearSystem:zeromatrix(ncoef, ncoef + ndepvar),
 146     for idat:1 thru ndat do (
 147       TransformedData:makelist(product(if PowersList[icoef][ivar] = 0 then 1
 148                                        else Mat[idat,ivar]^PowersList[icoef][ivar],
 149                                        ivar, 1, nvar),
 150                                icoef, 1, ncoef),
 151       for irow:1 thru ncoef do (
 152         for icol:1 thru ncoef do
 153           LinearSystem[irow,icol]:LinearSystem[irow,icol] +
 154                             TransformedData[irow] * TransformedData[icol],
 155         for ivar:1 thru ndepvar do
 156           LinearSystem[irow, ncoef + ivar]:LinearSystem[irow, ncoef + ivar] +
 157                                      DepCols[ivar][idat][1] * TransformedData[irow]
 158       )
 159     ),
 160     for ivar:1 thru ndepvar do
 161       DepSums[ivar]:col(LinearSystem, ncoef+ivar), /* save this info before modifying it */
 162
 163     /* Calculation of polynomial coefficients by solving the linear system with
 164        the Gauss method */
 165     PolyCoef:zeromatrix(ndepvar, ncoef),
 166     LinearSystem:ev(triangularize(LinearSystem), keepfloat:true),
 167     if product(LinearSystem[icoef,icoef], icoef, 1, ncoef) = 0 then (
 168       print("plsquares: insufficient number of independent data rows."),
 169       return(false)
 170     ),
 171     for ivar:1 thru ndepvar do (
 172       for irow:ncoef thru 1 step -1 do (
 173         PolyCoef[ivar,irow]:LinearSystem[irow,ncoef+ivar],
 174         for icol:irow+1 thru ncoef do
 175           PolyCoef[ivar,irow]:PolyCoef[ivar,irow] -
 176                                LinearSystem[irow,icol] * PolyCoef[ivar,icol],
 177         PolyCoef[ivar,irow]:PolyCoef[ivar,irow] / LinearSystem[irow,irow]
 178       )
 179     ),
 180
 181     /* Calculation and display of the determination coefficient(s) */
 182     DETCOEF:makelist(1 -
 183                      (sum(DepCols[ivar][idat][1]^2, idat, 1, ndat) -
 184                       sum(PolyCoef[ivar,icoef] * DepSums[ivar][icoef][1],icoef,1,ncoef)) /
 185                      (sum(DepCols[ivar][idat][1]^2, idat, 1, ndat) -
 186                       sum(DepCols[ivar][idat][1],idat,1,ndat)^2 / ndat),
 187                  ivar, 1, ndepvar),
 188     if atom(depvars) then DETCOEF:DETCOEF[1],
 189     print("      Determination Coefficient for", depvars, "=", float(DETCOEF)),
 190
 191     /* Construction and return of the polynomial(s) */
 192     PolyList:makelist(DepVarList[ivar]=
 193                       xthru(sum(PolyCoef[ivar,icoef] *
 194                                 product(VarList[jvar]^PowersList[icoef][jvar],
 195                                          jvar, 1 ,nvar),
 196                                 icoef, 1, ncoef)),
 197                       ivar, 1, ndepvar),
 198     if atom(depvars) then return(PolyList[1])
 199                      else return(PolyList)
 200   )$