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

[John Pye]

Hi Patrick,

Patrick Alken wrote:
> One solution to your problem is to not use the function
> gsl_multifit_linear, and instead call the gsl_linalg_svd_decomp routine
> yourself, search for singular values in the matrix S and when you find
> them the corresponding vector in V is a null space vector.
>   
Ok, I've used the gsl_linalg_SV_decomp function to obtain the singular
value decomposition, and indeed I have a very small (1e-33) value in my
vector 'S'.

I'm still feeling my way around GSL: what do you suggest I do next to
find the correct fitting plane? Can I use the remaining values in 'S' to
form my solution? Or do I need to set up another linear fit with a
reduced set of variables?

BTW if it were possible to fix gsl_multifit_linear for this case (or
alternatively to provide gsl_multifit_mul) that would be great!

Cheers
JP

> On Tue, Nov 20, 2007 at 01:07:03PM +1100, 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
>>> http://lists.gnu.org/mailman/listinfo/help-gsl
>>>       
>> _______________________________________________
>> Help-gsl mailing list
>> address@hidden
>> http://lists.gnu.org/mailman/listinfo/help-gsl
>>
>>     
>
>
> _______________________________________________
> Help-gsl mailing list
> address@hidden
> http://lists.gnu.org/mailman/listinfo/help-gsl
>