[Help-gsl] Nonlinear Curve Fitting aways returns "succes"

[Chris]

Hey everyone,
I'm having some issues with getting fitting a nonlinear model to some
data... When I call the gsl_multifit_fdfsolver_iterate() function, it ALWAYS
returns successful after the first iteration, no matter how far off it is.
I have a 5 parameter model and I believe I have set everything up correctly.
Here is some code for the main section that initializes the solver:

[code]

    //
    // Setup
    //
    int status = GSL_CONTINUE;
    int iterationCount = 0;
    int dataSize = (int)gsr->mRawGSRData.size();
    gsl_vector *initialGuessVector, *resultsVector;

    //
    // 1. Initialize solver using LMDER (Levenberg-Marquardt) for 'dataSize'
data points and a 5 parameter model...
    //
    const gsl_multifit_fdfsolver_type* T = gsl_multifit_fdfsolver_lmder;
    gsl_multifit_fdfsolver* fitSolver = gsl_multifit_fdfsolver_alloc (T,
dataSize, 5);

    //
    // 2. Create new multifit function data object
    //
    gsl_multifit_function_fdf fitFunction;

    // Model
    fitFunction.f = &(GSRMetricsCalculator::FittingFunction);

    // Jacobian
    fitFunction.df = &(GSRMetricsCalculator::JacobianMatrix);

    // Function to call f(x) and df(x)
    fitFunction.fdf = &(GSRMetricsCalculator::CalculateFitAndMatrix);

    // the number of functions, i.e. the number of components of the vector
f.
    fitFunction.n = dataSize;

    // the number of independent variables, i.e. the number of components of
the vector x.
    fitFunction.p = 5;

    // a pointer to the arbitrary parameters of the function.
    fitFunction.params = &gsr->mRawGSRData;

    //
    // 3. Create the initial guess vector
    //

    initialGuessVector = gsl_vector_alloc(5);
    gsl_vector_set(initialGuessVector, TIME_DECAY, 0.5);    // decay-time
constant
    // etc. etc..

    //
    // 4. Set the solver
    //
    gsl_multifit_fdfsolver_set(fitSolver, &fitFunction, initialGuessVector);

    //
    // 5. Run iterations
    //
    while((status == GSL_CONTINUE) && (iterationCount < 500))
    {
        status = gsl_multifit_fdfsolver_iterate(fitSolver); // always
returns "GSL_SUCCESS"
        iterationCount++;
    }

[/code]

Now, at step 5, no matter how far off I throw the initial guesses the
algorithm always seems to return "GSL_SUCCESS" (0) and never returns
"GSL_CONTINUE" (-2). Am I missing something here?? Below is the code for my
fitting and matrix functions (if needed):

FITTING FUNCTION

[code]

    //
    // Extract parameters from the vector "x" and set up variables.
    //

    double t_decay      = gsl_vector_get (x, TIME_DECAY             );  //
decay-time constant
    double t_rise       = gsl_vector_get (x, TIME_RISE              );  //
rise-time constant
    double gain1        = gsl_vector_get (x, GAIN_1                 );  //
gain for SCR 1
    double t_onset1     = gsl_vector_get (x, TIME_ONSET_1           );  //
onset-time for SCR 1
    double skinConLevel = gsl_vector_get (x, SKIN_CONDUCTANCE_LEVEL );  //
skin conductance level

    double SCR1;
    vector<double>* gsrData = (vector<double>*)(dataToFit);

    //
    // Populate the result vector, "f", using the 5-parameter model
described in:
    // "Decomposing skin conductance into tonic and phasic components"
    //

    int sizeOfData = (int)gsrData->size();
    for( int i = 0; i < sizeOfData; i++)
    {
        double t = (*gsrData)[i];

       // (calculations removed to conserve space in this email)

        gsl_vector_set (f, i, result);
    }

    return GSL_SUCCESS;

[/code]

MATRIX FUNCTION

[code]

    //
    // Store parameters
    //

    double t_decay      = gsl_vector_get (x, TIME_DECAY);               //
decay-time constant
    double t_rise       = gsl_vector_get (x, TIME_RISE);                //
rise-time constant
    double gain1        = gsl_vector_get (x, GAIN_1);                   //
gain for SCR 1
    double t_onset1     = gsl_vector_get (x, TIME_ONSET_1);             //
onset-time for SCR 1
    double skinConLevel = gsl_vector_get (x, SKIN_CONDUCTANCE_LEVEL);   //
skin conductance level

    vector<double>* gsrData = (vector<double>*)(dataToFit);

    //
    // Populate the jacobian matrix "j" by taking the partial derivative
with
    // respect to each parameter (columns) for every point of data (rows)
    //

    int sizeOfData = (int)gsrData->size();
    for( int i = 0; i < sizeOfData; i++)
    {
        double t = (*gsrData)[i];

        // (Calculations removed to conserve space in this email)

        // Store results in row "i" of the jacobian matrix
        gsl_matrix_set (j, i, TIME_DECAY            , d1);
        gsl_matrix_set (j, i, TIME_RISE             , d2);
        gsl_matrix_set (j, i, GAIN_1                , d3);
        gsl_matrix_set (j, i, TIME_ONSET_1          , d4);
        gsl_matrix_set (j, i, SKIN_CONDUCTANCE_LEVEL, 1);
    }

    return GSL_SUCCESS;

[/code]

I appreciate the help in advance.

~ Chris