Minpack added for solving non-linear equations and non-linear least-squares AND linear fitting

On 19. Jan, 2009, at 8:01 PM, Raymond Toy wrote:
>
> Good suggestion.  But optimization is such a vague name.  Does it load
> zillions of packages that are remotely related to optimization?  Even
> if the user only needs one simple optimization routine?
>
>     ziga> Hope someone writes nice wrapping functions for minpack.
>
> That someone could be you! :-)
>
> Ray

Maybe load(minimisation) is more appropriate.. the package would  
contain f2cl minpack and wrapper maxima level functions...

Yes, maybe I'll write it :) But I'm not yet proficient in lisp  
(learning it for a month or two, examining maxima code), however I've  
invested cca. 1 "man week" of work into wxMaxima interface since it's  
written in C++ which I know relatively well. So I've in some way  
contributed to Maxima, though not directly :)

While speaking of contributions, here is my user level, super easy to  
use, linear fitting function:

/******** MAXIMA CODE ************/
fit(data, functionlist, vars, [optional_args]) := block(
     [numer:true,float:true,
     lambdalist, design_matrix,
     alpha_matrix, beta_vector, sol,
     vector_b, out_format, out_list, _model, _out_subs:[] ],

/* PARAMETER: vars */
/* check if we have a single variable in vars, make it a list */
if atom(vars) then vars:[vars],

/* PARAMETER: data */
/* convert data to a matrix if it's a list */
if listp(data) then data: apply('matrix, data) ,
/* check length of datapoints against the number of vars */
if ( is(length(data[1]) < (length(vars) + 1)) or is(length(data[1]) >  
(length(vars) + 1))) then
     error("fit: The number of variables does not match data."),

/* PARAMETER: [optional_args] */
/* search [optional_args] for 'output = , default 'model */
out_format: assoc( 'output, optional_args, 'model),
/* out_list used for "member" search */
out_list: flatten([out_format]),

/* PARAMETER: functionlist */
/* Convert functionlist to a list of lambda functions with vars  
arguments */
lambdalist: map(
lambda([what], apply( 'lambda, what) ),
makelist( [vars,functionlist[function_index]] , function_index, 1,  
length(functionlist))
                 ),

/* FITTING */
/* Construct a design matrix (evaluate lambdalist on all datapoints) */
design_matrix:
     makelist(
         makelist(
             apply(lambdalist[function_index],
                   makelist(data[data_index,k], k, 1 , length(vars)) )
         , function_index, 1, length(functionlist))
     , data_index, 1, length(data))
,
design_matrix: apply('matrix, design_matrix) /* convert listoflist to  
matrix */
,
/* save design_matrix to _out_subs if user requested it in out_list */
if member('design_matrix, out_list) then
     (_out_subs:cons( 'design_matrix = design_matrix,  _out_subs)),

vector_b:
makelist(
data[data_index, length(vars) + 1]
, data_index, 1, length(data))
,
alpha_matrix: transpose(design_matrix) . design_matrix,
beta_vector: transpose(design_matrix) . vector_b,
sol:  alpha_matrix^^-1 . beta_vector,
sol: subst( '"[", 'matrix , transpose(sol) )[1],

/* OUTPUT FORMULATION */
/* write fit information to _out_subs */
if member('parameters, out_list) then
     (_out_subs:cons( 'parameters = sol , _out_subs)),
_model: apply("+" , sol * functionlist ),
if member('model, out_list) then
     (_out_subs:cons( 'model = _model , _out_subs)),
if member('residuals, out_list) then
     (_out_subs:cons( 'residuals = fit_residuals(data,_model,  
vars) ,_out_subs))
,
subst(_out_subs , out_format)
) $

fit_residuals( data, fitfunction, vars ) := block(
[numer:true,lambdafun],
lambdafun: apply( 'lambda, [vars , fitfunction]),
makelist(
data[data_index][length(vars) + 1] -
apply( lambdafun, makelist( data[data_index][k],k,1,length(vars) ))
, data_index, 1, length(data))
) $

/********* MAXIMA CODE ***********/

Example of use:

/********* MAXIMA CODE ***********/
dejta: [ [1.0, 2.0], [ 2.0, 3.0], [4.0, 5.0], [8.0,8.0], [6.0,5.0],  
[12.0,-5]];
/* basic usage */
fit(dejta, [1, x], x);

/* 'expert' usage */
fit(dejta, [1, x], [x], output = '[model,residuals, 
[design_matrix,parameters]]);

/* 2d scenario */
data2d: [ [1, 2, 2] , [ 3, 2, 5], [1, 5, 4], [3, 5, 5],[2,1,3], 
[-1,0,7]] $
fit2d: fit( data2d, [1,x,y,x*y], [x,y] );

/********* MAXIMA CODE ***********/

Function fit accepts either a list of lists or a matrix. By default  
it returns just the model, but the user can specify it's output,  
since there is a lot of information one can extract from a fit. Right  
now "fit" can output only
model, residuals, parameters, design_matrix
but the function can be extended to calculate various information  
about the fit. I'm not a fitting expert, so I'm open to suggestions  
which statisics can be extracted from a fit (anyone?). I will add  
"chi_squared" output element which is the sum of residuals squared  
(or is that called R^2 ?) and support for sigmas (error estimates of  
the data).
I intend to write a nonlinear fitting function with similar calling  
sequence with support for multiple datasets (I've read that someone  
mentioned it on the mailing list).

I think configurable output of a fitting function is a nice touch  
maxima can provide over MATLAB and such. It makes "fit" much more  
powerful, while retaining it's ease of use, when output is not  
specified. User might not be concerned about the model itself but  
some statistic. One can compare the quality of the fit for different  
polynomials easily:

somedata: [[1,1],[2,2],[3,-1],[5,5],[-1,6],[8,8]]$ /*something*/
makelist([imax, apply("+",   fit( somedata, makelist(x^i, i, 0,  
imax) , x, 'output = 'residuals)^2 ) ], imax, 1, 3);

I think such configurable output really shows the power of Maxima  
compared to systems like MATLAB. Maxima numerical functions should  
make maximum use of it's symbolic core to gain an edge over faster  
but not as flexible numerical systems. Lisp vs Fortran is similar in  
some ways to Maxima/CAS vs MATLAB.

I'm open for suggestions about possible outputs of a fitting  
function... (preferably with some link to resources regarding  
implementation/ calculation of these outputs).

Ziga