Re: [Help-gsl] puzzled with result from gsl_multifit_linear...

[John Gehman]

Hi John,

Perhaps you're getting sucked in to your own simplified example in  
trying to get some code working? Presumably the bigger picture is  
that you have many more than three points, each with some error  
associated with them, and you're trying to fit the best plane through  
them. In this case your objective function to be minimized through  
optimisation of a,b,c would be something like the square of the  
distances between a given trial plane and each of your points, akin  
to a line or other curve in two dimensions. I don't know if the  
linear regression functions are extendable to three-space, but I  
would guess you could easily use the nonlinear fitting functions.

Otherwise, if this is more a geometry problem than a model-fitting  
problem, while it' s been some time since I thought about this, I  
somewhat remember settling on easily defining two vectors from three  
points, and computing the cross product which will be the normal to  
the plane, and the coefficients of that normal vector bears some  
clear relationship to the a,b,c,d of the plane (obviously I don't  
recall very well, but should be easy to find).

Maybe this helps?

Cheers,
john



On 20/11/2007, at 2:38 PM, John Pye wrote:

> G'day to you too...
>
> Thanks for the reply, John. I think you're right. In general a plane
> requires the 'd' parameter, but only really for the degenerate case
> where the plane passes through the origin. For any other case, the
> equation can be rearranged such that the 'd' value is replaced by a  
> '1'
> on the RHS.
>
> This lead me to wonder how I can implement a plane-fitting routine  
> that
> is robust and general. I think that the correct approach might be to
> offset all the points by the centroid, and then fit the equation
> ax+by+cz=0. There is a theorem about the centroid always lying on  
> the plane.
>
> So now I tried to implement this using GSL but it seems that GSL  
> lacks a
> function like 'gsl_multifit_mul'. When I try to fit my plane equation
> using gsl_multifit_linear, with y[i]=0 for all i, I get the  
> predictable
> result that a=0, b=0, c=0. So I don't think GSL can do what I need  
> it to
> do. Or perhaps there is a cunning workaround?
>
> Thoughts?
>
> Cheers
> JP
>
> John Gehman wrote:
> > G'day John,
> >
> > I haven't looked at the examples all that closely, but just  
> > quickly, I
> > believe a general equation for a plane is ax+by+cz+d=0, i.e. in your
> > equation you're effectually fixing d at -1, while I think it  
> > should be
> > d=0 (zero) in your example?
> >
> > Cheers,
> > john
> >
> > ---------------------------------------------------------
> >
> >
> > John Gehman
> >
> > Research Fellow
> >
> > School of Chemistry
> >
> > University of Melbourne
> >
> > Australia
> >
> >
> > On 20/11/2007, at 12:53 PM, John Pye wrote:
> >
> > > Hi all
> > >
> > > I have a question about the use of the gsl_multifit_linear  
> > > routine that
> > > perhaps is a question more about geometry/algebra than coding,  
> > > but I'm
> > > not sure.
> > >
> > > I want to construct a routine that fits a plane (in three  
> > > dimensions)
> > > through a set of data points (x,y,z). I have set up  
> > > gsl_multifit_linear
> > > to fit the plane equation a*x + b*y + c*z = 1to my data, and for  
> > > most
> > > cases that seems to work OK. However there are a few degenerate  
> > > cases
> > > that don't work, and I'm trying to work out what I should do. Is  
> > > there a
> > > better equation that describes a plane?
> > >
> > > The following example shows what happens when I try to fit a  
> > > plane to
> > > the points (0,0,0), (1,0,0), (1,1,1), (0,1,1). These points  
> > > represent an
> > > inclined rectangle with one side on the positive x-axis.
> > >
> > > I would expect the fit to these points to be 0*x + INF*y - INF*z  
> > > = 1, or
> > > else for it to give an error. But instead, gsl_multifit_linear  
> > > gives me
> > > the result 0.666*x + 0.333*y + 0.333*z = 1, which is wrong.
> > >
> > > I am obviously making a logic error here, but I'm a bit stuck on  
> > > what it
> > > is, and what the correct general way of working around it would  
> > > be? Can
> > > anyone here offer any suggestions?
> > >
> > > Cheers
> > > JP
> > >
> > >
> > > // for the fitting of a plane through a group of points, testfit.cpp
> > > #include <gsl/gsl_vector.h>
> > > #include <gsl/gsl_matrix.h>
> > > #include <gsl/gsl_multifit.h>
> > > #include <iostream>
> > >
> > > using namespace std;
> > >
> > > int main(void){
> > >
> > >     int n = 4;
> > >     int num_params = 3;
> > >
> > >     gsl_multifit_linear_workspace *work =  
> > > gsl_multifit_linear_alloc(n,
> > > num_params);
> > >
> > >     gsl_vector *y = gsl_vector_calloc(n);
> > >     for(unsigned i=0; i < n; ++i){
> > >         gsl_vector_set(y,i,1.0);
> > >     }
> > >
> > >     gsl_matrix *X = gsl_matrix_calloc(n, num_params);
> > >     gsl_vector *c = gsl_vector_calloc(num_params);
> > >     gsl_matrix *cov = gsl_matrix_calloc(num_params, num_params);
> > >
> > > #define SET(I,A,B,C) \
> > >     gsl_matrix_set(X,I,0,A); gsl_matrix_set(X,I,1,B);
> > > gsl_matrix_set(X,I,2,C)
> > >
> > >     SET(0, 0,0,0);
> > >     SET(1, 1,0,0);
> > >     SET(2, 1,1,1);
> > >     SET(3, 0,1,1);
> > >
> > >     double chisq;
> > >     int res = gsl_multifit_linear(X,y,c,cov,&chisq,work);
> > >
> > >     cerr << "FIT FOUND:" << endl;
> > >     for(int i=0;i<3;++i){
> > >         cerr << "  c[" <<i << "]=" <<gsl_vector_get(c,i) << endl;
> > >     }
> > >
> > > }
> > >
> > >
> > > _______________________________________________
> > > Help-gsl mailing list
> > > address@hidden <mailto:address@hidden>
> > > http://lists.gnu.org/mailman/listinfo/help-gsl