# [Maxima] Minpack added for solving non-linear equations and non-linear least-squares AND linear fitting

Žiga Lenarčič ziga.lenarcic at gmail.com
Mon Jan 19 14:05:24 CST 2009

```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

```