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

[John Pye]

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;
    }

}