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
http://lists.gnu.org/mailman/listinfo/help-gsl