Walking Randomly demonstrated in his blog how to solve non-linear fit problems with various software systems. I thought I might try to solve this in Euler.

We have the following data.

>xdata = [-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9];
>ydata = [0.699369,0.700462,0.695354,1.03905,1.97389,2.41143, ...
>  1.91091,0.919576,-0.730975,-1.42001];

and the non-linear model

\(y = p_1 \cos(p_2 x) + p_2 \sin(p_1 x).\)

We need to determine p so that the summed squares of the errors is minimal. Let us use Powell’s minimizer in EMT. This algorithm searches along several directions one after the other and then adjusts those directions. For that, it takes the current total descent into the set of directions. It is not a quick algorithm and cannot be used in high dimensional problems.

We define a function q(p) and a function qerr(p) for the error.

>function q(x,p) := p[1]*cos(p[2]*x)+p[2]*sin(p[1]*x);
>function qerr(p) := sqrt(sum((ydata-q(xdata,p))^2));
>p=powellmin("qerr",[1,0.2])
 [1.88185,  0.70023]
>plot2d(xdata,ydata,>points);
>plot2d("q(x,p)",>add,color=red):

That was easy. We could just as well have taken the Nelder-Mead algorithm. That algorithm adjusts a simplex of points so that the simplex walks towards the minimum. It is usually a good and stable algorithm.

>neldermin("qerr",[1,0.2])
 [1.88185,  0.70023]

It is interesting to use an algorithm with a derivative. The Levenberg-Marquardt method is such an algorithm. Since it minimizes functions from n dimensions to m dimensions with m>n, it walks in the space of n variables to minimize the norm of the function, which maps to m variables. For this, it minimizes the tangential function.

\(T(\tilde x) = f(x) + Df(x) (x-\tilde x)\)

The implementation in EMT needs the functions f(v) and Df(v).

>function f(p1,p2,x,y) &= y-p1*cos(p2*x)-p2*sin(p1*x)
                    y - p1 cos(p2 x) - p2 sin(p1 x)
>function fv(p) := f(p[1],p[2],xdata',ydata')';
>function df1(p1,p2,x) &= diff(f(p1,p2,x,y),p1);
>function df2(p1,p2,x) &= diff(f(p1,p2,x,y),p2);
>function dfv(p) := df1(p[1],p[2],xdata')|df2(p[1],p[2],xdata');
>nlfit("fv","dfv",[1,0.2])
 [1.88185,  0.70023]

The construction of dfv(p) for a vector x is a bit tedious, since EMT cannot combine a matrix from columns using the simple vector format [df1,df2] which the gradient function returns. This will be added in the next versions of EMT. Then the solution shortens as follows.

>function f(p1,p2,x,y) &= y-p1*cos(p2*x)-p2*sin(p1*x)

                    y - p1 cos(p2 x) - p2 sin(p1 x)

>function fv(p) := f(p[1],p[2],xdata',ydata')';
>function df(p1,p2,x) &= gradient(f(p1,p2,x,y),[p1,p2])

        [- cos(p2 x) - p2 x cos(p1 x), p1 x sin(p2 x) - sin(p1 x)]

>function dfv(p) := df(p[1],p[2],xdata');
>nlfit("fv","dfv",[1,0.2])
 [1.88185,  0.70023]

If we want to make the process easier for the user, we could try to write a function modelfit(), which accepts a model, a starting guess and the data. We need a function, which evaluates the error of the current parameter. The model function is passed by name in f$. We call neldermin() for this function to find the minimum. The semicolon parameters f$, x and y are passed from meldermin() to %evalmodel().

>function %evalmodel (p,f$,x,y) := sum((y'-f$(x',p))'^2);
>function modelfit (f$,p,x,y) := neldermin("%evalmodel",p;f$,x,y);

Now the user can simply call the model in the following way.

>function model(x,p) := p[1]*cos(p[2]*x)+p[2]*sin(p[1]*x);
>modelfit("model",[1,0.2],xdata,ydata)
 [1.88185,  0.70023]