lsquares

Hi,

I just found the lsquares package in maxima. I have couple of questions:

1) Is it possible to simultaneously fit multiple data sets with each 
having its own fitting function in such a way that they share parameters?
For example, I run into this sort of problem often in chemical kinetics 
where multiple components are followed experimentally and each have 
their own kinetic equation. The rate constants are shared between the 
two equations when fitting the data.

2) Are there routines to calculate r^2, covariance matrix, standard 
error estimates ? Of course, for non-linear least squares things are a 
bit tricky but many things can be done at least approximately.

I have been using the following maxima program for (non-linear) least 
squares. I would like to get rid of this and use the built-in routines 
in maxima.

Thanks,

Jussi Eloranta
Department of chemistry and biochemistry
Cal State Northridge

---

/*
 * Unconstrained nonlinear least squares fit.
 *
 * nlsq(funcs, ydata, xdata, yvars, xvars, params, initvars, errtol, debug)
 *
 * funcs    = list of functions for fitting [fun1, fun2, ...].
 * ydata    = list of y-data vectors [yvec1, yvec2, ...].
 * xdata    = list of x-data vectors [xvec1, xvec2, ...].
 * yvars    = list of y variables for the functions given in funcs [y1, 
y2, ...].
 * xvars    = list of x variables for the functions given in funcs [x1, 
x2, ...].
 * params   = list of parameters to be fitted [a, b, ...].
 * initvars = list of initial values for params [ai, bi, ...].
 * errtol   = requested error tolerance in BFGS (1E-4).
 * debug    = optimization output vector for BFGS ([1, 2] = full output, 
[-1 0] = no output).
 *
 * Returns: [F, R2, V, C, E].
 *
 * F  = list of optimized functions.
 * R2 = list of r^2 values for each function.
 * V  = list of optimized parameter values.
 * C  = covariance matrix for the parameters.
 * E  = list of standard errors for parameters.
 *
 */

load("lbfgs");

nlsq([ArgList]):=
  block([narg, funcs, ydata, xdata, yvars, xvars, params, initvars, errtol,
         debug, neqs, i, sq, tmp, lbfgs_nfeval_max:1000, /* BFGS tends 
to converge somewhat slowly... */
         fun, fit, F, ss, sstot, sm, st, R2, C, lp, E, sigma2, Vm, k, np],

  narg:length(ArgList),
  if narg # 9 then (
    print("nlsq: bad number of function arguments (requires 9 arguments)."),
    return(false)
  ),

  funcs:ArgList[1],
  ydata:ArgList[2],
  xdata:ArgList[3],
  yvars:ArgList[4],
  xvars:ArgList[5],
  params:ArgList[6],
  initvars:ArgList[7],
  errtol:ArgList[8],
  debug:ArgList[9],

  neqs:length(funcs),
  lp:length(params),

  if neqs # length(ydata) or neqs # length(xdata) or neqs # 
length(yvars) or neqs # length(xvars) then (
    print("nlsq: Inconsistent lengths for arguments 1 - 5."),
    return(false)
  ),

  /* Construct the least squares sum and call lbfgs to optimize */
  tmp:0,
  for i:1 thru neqs do (
    if length(xdata[i]) # length(ydata[i]) then (
      print("nlsq: Inconsistent lengths for X and Y in set ", i),
      return(false)
    ),
    sq:(lhs(funcs[i]) - rhs(funcs[i]))^2,
    sq:subst('xdata[i][j], xvars[i], sq),
    sq:subst('ydata[i][j], yvars[i], sq),
    tmp:tmp + 'sum(sq, j, 1, length(xdata[i]))
  ),
/*  fun:ev(tmp,nouns),*/
  fun:tmp,
  fit:lbfgs(fun, params, initvars, errtol, debug),
  if fit = [] then (
    print("nlsq: BFGS convergence error. Fit failed."),
    return(false)
  ),

  /* Substitute the optimized values back into the equations */
  F:makelist(0,i,1,neqs),
  for i:1 thru neqs do (
    F[i]:float(rhs(funcs[i])),
    for j:1 thru lp do F[i]:subst(rhs(fit[j]), params[j], F[i])
  ),

  /* Calculate r^2 for each data set */
  sstot:0,
  R2:makelist(0,i,1,neqs),
  for i:1 thru neqs do (
    ss:sum((subst(xdata[i][j], xvars[i], F[i]) - ydata[i][j])^2, j, 1, 
length(xdata[i])),
    sstot:sstot + ss,
    sm:sum(ydata[i][j], j, 1, length(xdata[i])) / length(xdata[i]),
    st:sum((subst(xdata[i][j], xvars[i], F[i]) - sm)^2, j, 1, 
length(xdata[i])),
    R2[i]:float(ev(1 - ss / st,nouns))
  ),

  /* Calculate the covariance matrix */
  C:ematrix(lp, lp, 0, 1, 1),
  for i:1 thru lp do
    for j:1 thru lp do (
      C[i,j]:diff(fun,params[i],1,params[j],1),
      for k:1 thru lp do C[i,j]:subst(rhs(fit[k]),params[k],C[i,j])
    ),
  C:float(ev(C,nouns)),
  /* watch out when lp = 1, invert returns a number not a matrix */
  if lp = 1 then (
    C[1,1]:1/C[1,1]
  ) else (
    C:invert(C)
  ),

  /* Calculate standard errors */
  np:sum(length(xdata[i]), i, 1, neqs),
  sigma2:sstot / (np - lp),
  E:makelist(0,i,1,lp),
  V:makelist(0,i,1,lp),
  for i:1 thru lp do (
    E[i]:float(sqrt(sigma2 * C[i,i])),
    V[i]:float(rhs(fit[i]))
  ),

  return([F, R2, V, C, E])
);