# 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
> restart:
> read(`/home/mat331/lsq_data`);
defined line_pts(), quadratic_pts(), cubic_pts(), and circle_pts()
> cubdata := cubic_pts():
> 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((cub(cubdata[i][1],a,b,c,d) - cubdata[i][2]
> )^2,                        i=1..nops(cubdata));

G := (a, b, c, d) ->

    nops(cubdata)
        -----
         \                                                        2
          )       (cub(cubdata[i][1], a, b, c, d) - cubdata[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
    -.3688032502 x  + .3920904491 x  + 2.656274173 x + 1.999084485

# Finally, lets make a picture, just to see how well it fits
> with(plots):
> cplot:= plot(cub(x,a,b,c,d), x=-3..3, 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:
> 
> 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 = -.3688032508 x  + .3920904505 x  + 2.656274177 x + 1.999084480

# they match to 8 figures, which is all we could have hoped for.  But of
# course, they should.