Fitting in Maxima



Maxima is actually lacking a proper fitting facilities (like many  
other things). lsquares is a poor excuse for a fitting package in  
many aspects:
-slow
-complicated to use
-unintuitive
-cryptic name (should be LeastSquares/least_squares, not lsquares_mse  
and what not)
I think the latter three (if not all four) aspects of lsquares  
package applies well to a lot of maxima functions, if they even exist  
(still no ODE integrator or fast numerical linear algebra  
functions)... Especialy cryptic names are the plague of maxima. What  
can a user dechipher from 'create_list' and 'makelist'. It makes no  
sense..

However not to be too critical towards maxima developers (since a lot  
of problems really originate from 30 years old code, which nobody  
dares to rewrite), I find the new draw package extremely well  
designed (esp. compared to plot2d), easy to use and useful (though  
some names, like "terminal" or "enhanced3d", originating from gnuplot  
should be translated to meaningful names like "output".) and I really  
like the fact that new gamma implementations feel like a complete  
coherent package, naming and functionality-wise.

There sadly is no "manual" written for extending maxima, no  
directions for adding functionality to maxima. I feel if new  
additions to maxima aren't guided by some central guidance document,  
different maxima level naming conventions, internal naming  
conventions, packages doubling or only partially implementing desired  
featureset ('vect' and 'vector' for instance) will keep maxima being  
a nonelegant mess of poor implementations (which you must admit in  
many areas it is).
I would really like, somebody wrote such a document with maybe some  
general plan of maxima features to be added or reimplemented in a  
better way.

Of course I'm not only complaining but I will gladly contribute. I  
have very little experience with lisp and symbolic computations  
source code, but have some experience with numerical methods..

I will try to add linear fitting and also nonlinear fitting  
functionality to maxima via a nice full featured package. I'm also  
thinking of writing numerical ode solving package, but that's a much  
tougher task, so maybe in the future.

I have a working linear fitting function written in Maxima language  
and while writing it, many questions regarding adding functionality  
to Maxima.

1) What is the current naming convenction for new functions in Maxima?
It seems that new functions are named in a  
whole_words_separated_with_underscore. I personally much prefer the  
Mathematica's CapitalsForEachWord, which is easier to read, easier to  
type and also leaves more room for user defined variables or  
functions (which are typically small caps). I understand the need for  
backwards compatibility, but new conventions could be applied to new  
functionality easily. I'm glad that at least MATLAB's abrvtdfunnames  
aren't encouraged, but are sadly very common in maxima.

2) Are there any internal (lisp) naming conventions? Code reuse would  
benefit from this.

3) When should new functionality be implemented in Maxima language  
and when in lisp if for instance some problem can be solved in both?

4) Will there be possibility to call compiled numerical libraries  
(ATLAS ..) from lisp in the future? I know gcl is kinda the limiting  
factor right now..  I think f2cl conversions of numerical libraries  
should not be a part of maxima. Especially if the user has to call  
them with their fortran names. It's better to implement an  
unoptimised original algorithm in lisp by hand (and tailored to  
maxima's needs), since f2cl translations do not inherit the speed of  
original fortran porgrams, neither are they consistent with other  
maxima code or gain anything. So it's like using fortran without  
speed, why would you do that? Ideally we should use ATLAS for  
numerical linear algebra, so speed-wise maxima would be on par with  
MATLAB and mathematica. Any chance of this becoming reality?

5) List of lists vs matrix : transpose does not work on list of  
lists, therefore one has to convert list of lists to matrix,  
transpose and convert back. I don't see what functionality is gained  
by reprezenting matrices with a 'matrix. In case we use a special  
construct for matrices, we should at least make something of it. A  
matrix could carry some useful flags about it's contents: real/ 
complex, diagonal, symetric, numerical(all elements) ... which would  
then updated at any operation on matrix and used when solution or a  
system or eigenvalues are needed for choosing the best algorithm for  
the job (i.e. if all elements are numeric, the matrix is passed to  
ATLAS numerical libraries). Transpose should accept also a list of  
lists. I'd rather see, 'matrix didn't exist since it's more trouble  
for the user, than what it's worth.

6) There is no real SVD in Maxima? (don't suggest lapack dgesvd)  
Maybe it should be implemented..

6) What functionality should Maxima provide regarding function fitting?
I intend to write functions with similar calling structure as  
Mathematicas
Fit - linear fitting (via missing SVD)
FindFit - nonlinear fitting via levenberg-marq. minimisation (from  
Numerical Recipes)
Suggestions regarding implementation of fitting are welcome!


Here is a linear fitting function in maxima language (hey! it works!)  
and examples:
/****** MAXIMA *******/
Fit(data, functionlist, vars) := block(
     [numer:true,float:true,
     lambdalist, design_matrix,
     alpha_matrix, beta_vector, sol,
     vector_b],

lambdalist: map(
lambda([what], apply( 'lambda, what) ),
makelist( [vars,functionlist[function_index]] , function_index, 1,  
length(functionlist)) )
,
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))
,
vector_b:
makelist(
data[data_index][length(vars) + 1]
, data_index, 1, length(data))
,
design_matrix: apply('matrix, design_matrix),
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]
,
apply("+" , sol * functionlist )
) $

/* 2d example */
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] );
load(draw) $
draw3d(
color = red,
point_size = 3, point_type = 7,
points(data2d),
color = black,
explicit( fit2d , x, -1, 4, y, 0, 5)
)$

/* 1d example, 1000 points */
dejta: makelist( [float(i/50), sin(i/50)+random(0.2)] , i, 0, 1000),  
numer$
fit:Fit( dejta, [1, x, sin(x)], [x]);
draw2d(
color = blue,
points(dejta),
color = red,
explicit( fit, x, dejta[1][1]-1, dejta[length(dejta)][1]+1)
) $


/******** MAXIMA ********/

I would like Fit to use svd instead of ^^-1, since it's more robust,  
but looks like i have to implement svd first? How and where should  
svd be implemented in Maxima?


hny
?iga Lenar?i?