Fitting a cubic to data

Let's now try to fit a cubic polynomial to some data. Again, maple's built-in package can do this for us, but we like to work harder than we need to, so we understand more. First, let's get some data, using the provided function

cubdata := read(`/home/mat331/lsq_data`):
plot(cubdata,style=point, axes=boxed);

As in the linear case, we will do the work "by hand".

To simplify typing, let's define the type of function we want to fit. We know (somehow) that our goal is to fit a cubic polynomial, so

 cub := (x,a,b,c,d) ->  a*x^3 + b*x^2 + c*x + d;

                                           3      2
              cub := (x, a, b, c, d) -> a x  + b x  + c x + d

Now we define an "error functional" which gives us the distance from our cubic to our data. Note that we use maple's "nops" command so that we don't need to explicitly say how many points there are.

G := (a,b,c,d) -> sum((cubdata(pts[i][1],a,b,c,d) - cubdata[i][2] )^2, 
                    i=1..nops(cubdata));

                    nops(cubdata)
                      -----
                       \                                                    2
  G := (a, b, c) ->     )     (cub(cubdata[i][1], a, b, c,d) - cubdat[i][2])
                       /
                      -----
                      i = 1

Now we just ask maple to find the (unique) critical point of G, which is our minimum. Since we used assign to set the coefficients, we won't see a result.

assign(solve( {diff( G(a,b,c,d), a)=0,
               diff( G(a,b,c,d), b)=0,
               diff( G(a,b,c,d), c)=0},
               diff( G(a,b,c,d), d)=0},
             {a,b,c,d} ));

but we can readily see what the answer is:

cub(x,a,b,c,d);

           3                2
          x  - .4806519894 x  - 10.94459406 x + 10.28284644

Finally, lets make a picture, just to see how well it fits

with(plots):
cplot := plot(cub(x,a,b,c,d), x=-4..4, axes=boxed):
pplot:= plot(cubdata, style=point):
display({cplot, pplot});

Just for comparison, we can let maple do this for us, using the built-in least-squares package:

 with(stats):
 fit[leastsquare[[x,y],y=A*x^3+B*x^2+C*x+D]]
        ([[cubdata[i][1]$i=1..nops(cubdata)],
          [cubdata[i][2] $i=1..nops(cubdata)]]);

              3                2
      y =    x  - .4806519892 x  - 10.94459410 x + 10.28284642