lsquares fitting and constrains for parameters (looking for positive values of parameters)

> If you do not need a symbolic solution,

You're right. I need concrete values. (I've forgot to say it in
problem statement).

> you might be interested in
> fmin_cobyla, which will produce the minimum of a function subject to
> equality or inequality constraints.

fmin_cobyla, I haven't knew about it. (btw. I see it's implemented in
Scipy as well.)

I try to figure out how to make it to obey constraints, with no
success. Results are not positive, and even worse than
lsquares_estimate.

As I doing sth wrong ?
Or sum equations are not fmin_cobyla's purpose, are they ?

To show clearly what is my problem, I've prepared this toy example ,
benchmarking lsquares with fmin_cobyla:
____________________________
display2d:false$
load(lsquares)$
load(fmin_cobyla)$
DataM: matrix( [1,2], [2,4], [3,6], [4,8] )$

D_vars: [x,y]$
pattern: y = a3*x**3 + a2*x**2 + a1*x + a0$
parameters: [a3,a2,a1,a0]$
mse : lsquares_mse (DataM, D_vars, pattern);

lsq_val: lsquares_estimates (DataM, D_vars, pattern, parameters);
cob_val: fmin_cobyla(mse, parameters, [1,1,1,1], constraints = [a3>=0,
a2>=0, a1>=0, a0>=0], iprint=1);

lsq_expr: ev(pattern,lsq_val);
cob_expr: ev(pattern,cob_val[1]);

quality_of_lsq: ev(mse,lsq_val),nouns;
quality_of_cob: ev(mse,cob_val[1]),nouns;
____________________________

Here is how it works:
____________________________
$ maxima -b /tmp/t.mac 2>/dev/null
Maxima 5.24.0 http://maxima.sourceforge.net
using Lisp SBCL 1.0.51
Distributed under the GNU Public License. See the file COPYING.
Dedicated to the memory of William Schelter.
The function bug_report() provides bug reporting information.
(%i1)                          batch(/tmp/t.mac)
(%i2)                          display2d : false
(%i3) load(lsquares)
(%i4) load(fmin_cobyla)
(%i5) DataM:matrix([1,2],[2,4],[3,6],[4,8])
(%i6) D_vars:[x,y]
(%i7) pattern:y = a0+a1*x+a2*x^2+a3*x^3
(%i8) parameters:[a3,a2,a1,a0]
(%i9) mse:lsquares_mse(DataM,D_vars,pattern)
(%o9) ('sum(('DataM[i,2]-a3*'DataM[i,1]^3-a2*'DataM[i,1]^2-a1*'DataM[i,1]-a0)
             ^2,i,1,4))
       /4
(%i10) lsq_val:lsquares_estimates(DataM,D_vars,pattern,parameters)
(%o10) [[a3 = 0,a2 = 0,a1 = 2,a0 = 0]]
(%i11) cob_val:fmin_cobyla(mse,parameters,[1,1,1,1],
                           constraints = [a3 >= 0,a2 >= 0,a1 >= 0,a0 >= 0],
                           iprint = 1)

   Return from subroutine COBYLA because the MAXFUN limit has been reached.

   NFVALS = 1000   F = 3.415514E-02    MAXCV = 5.648753E-20
   X =-5.648753E-20   1.754386E-01   1.073534E+00   1.018859E+00
(%o11) [[a3 = -5.648753124436182e-20,a2 = .1754386264111758,
         a1 = 1.073533597507375,a0 = 1.018858888779015],0.0341551363982473,
        1000,1]
(%i12) lsq_expr:ev(pattern,lsq_val)
(%o12) y = 2*x
(%i13) cob_expr:ev(pattern,cob_val[1])
(%o13) y = -5.648753124436182e-20*x^3+.1754386264111758*x^2
                                     +1.073533597507375*x+1.018858888779015
(%i14) ev(quality_of_lsq:ev(mse,lsq_val),nouns)
(%o14) 0
(%i15) ev(quality_of_cob:ev(mse,cob_val[1]),nouns)
(%o15) 0.0341551363982473
____________________________


Greg