[Help-gsl] Fit data values to gumbel distribution with GSL Levenberg-Marquardt

[Jochen]

Hi list,

after I found out, that it is quite hard (for me) to get gsl-curve-fitting done 
and I am not the only one, I created one example I want to push step by step.
Thanks in advance for comments which guide me in the right direction. As you 
can see, I got 30 values, I want 'mue' and 'beta' for the cumulated gumbel
distribution.  

Okay my code-skeleton follows, I hope it's correct and now I will try to fill 
the gaps. My main references are:
Doc: gnu.org/software/gsl/manual/html_node/index.html#Top
Gumbel-Distribution: en.wikipedia.org/wiki/Gumbel_distribution
Threads: lists.gnu.org/archive/html/help-gsl/2007-12/msg00004.html; 
                lists.gnu.org/archive/html/help-gsl/2005-12/msg00049.html

#include <fstream>
#include <gsl/gsl_histogram.h>
#include <gsl/gsl_multifit_nlin.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_sf_bessel.h>
#include <gsl/gsl_statistics.h>
#include <gsl/gsl_vector.h>
#include <math.h>
#include <iostream>
#include <stdlib.h>  
#include <stdio.h>
#include <time.h>
#include <vector>
//#include "nr3.h"
//#include "sort.h"

using namespace std;

struct curvefitdata 
{
    size_t n;
    double * y;
    double * mue;
    double * beta;
};

void curvefit_f(){}; // TODO
void curvefit_df(){}; // TODO
void curvefit_fdf(){}; // TODO
void curveFitPrintState(size_t iter, gsl_multifit_fdfsolver * s) {
    printf("iter: %3u x = % 15.8f % 15.8f % 15.8f |f(x)| = %g\n", iter, 
gsl_vector_get(s->x, 0), gsl_vector_get(s->x, 1), gsl_vector_get(s->x, 2), 
gsl_blas_dnrm2(s->f));
}

int main(int argc, char** argv) 
{
    // IDE-Properties:
    // Netbeans-Linker-Additional-Options: -lgsl -lgslcblas -lm
    // inclusions (hope I mentioned all here): 
        // #include <gsl/gsl_histogram.h>
        // #include <gsl/gsl_multifit_nlin.h>
        // #include <gsl/gsl_randist.h>
        // #include <gsl/gsl_rng.h>
        // #include <gsl/gsl_vector.h>
        
    cout << "GSL Curve Fitting-Example \n" << endl;
    
    const gsl_rng_type * T; // create NumberGenerator
    gsl_rng * r;
    gsl_rng_env_setup();
    T = gsl_rng_default;
    r = gsl_rng_alloc(T);
    
    // data-vector, 30 vals, maybe gumbel-distributed: 
    vector<double> uvalsRealVector;
    uvalsRealVector.push_back(0.0603461);
    uvalsRealVector.push_back(0.1436600);
    uvalsRealVector.push_back(0.1091960);
    uvalsRealVector.push_back(0.1441190);
    uvalsRealVector.push_back(0.1383480);
    uvalsRealVector.push_back(0.0634344);
    uvalsRealVector.push_back(0.0964405);
    uvalsRealVector.push_back(0.1011820);
    uvalsRealVector.push_back(0.1142970);
    uvalsRealVector.push_back(0.1221080);
    uvalsRealVector.push_back(0.1118600);
    uvalsRealVector.push_back(0.0852936);
    uvalsRealVector.push_back(0.1518220);
    uvalsRealVector.push_back(0.0780212);
    uvalsRealVector.push_back(0.0929122);
    uvalsRealVector.push_back(0.1219700);
    uvalsRealVector.push_back(0.1128130);
    uvalsRealVector.push_back(0.1478180);
    uvalsRealVector.push_back(0.0774807);
    uvalsRealVector.push_back(0.0975139);
    uvalsRealVector.push_back(0.0904746);
    uvalsRealVector.push_back(0.0898072);
    uvalsRealVector.push_back(0.1010630);
    uvalsRealVector.push_back(0.0739358);
    uvalsRealVector.push_back(0.0525509);
    uvalsRealVector.push_back(0.1371950);
    uvalsRealVector.push_back(0.1028160);
    uvalsRealVector.push_back(0.1035900);
    uvalsRealVector.push_back(0.1361230);
    uvalsRealVector.push_back(0.0845199);
    
    gsl_vector * uvalsReal = gsl_vector_alloc(30); // change data to 
GSL-expected data type
    for( int i = 0; i < 30; i++ )
    {
        gsl_vector_set( uvalsReal, i, uvalsRealVector.at( i ) );
    }
    
    // Curvefitting to Gumbel-Dist, Levenberg-Marquardt
    int status;         // solver-status
    int i;              // index
    int iter = 0;
    int constn = 40;    // number of observations, 'N' in Documentation
    const size_t n = constn;    
    const size_t p = 2;         // parameters
    const gsl_multifit_fdfsolver_type * T = gsl_multifit_fdfsolver_lmder;       
// solver-creation part 1
    gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc(T, 30, p);        
// solver-creation part 2
    printf("s is a '%s' solver\n", gsl_multifit_fdfsolver_name(s));             
// point to solver, solver-name
    gsl_multifit_function_fdf f;        // function
    double x_init[3] = {0.2, 0.002};    // starting guesses
    gsl_matrix *covar = gsl_matrix_alloc(p, p); // matrix for paras
    
    // data types, functions for curve-fitting-procedure (struct, f, df, fdf)
    double y[constn];
    double mue[constn];
    double beta[constn];
    struct curvefitdata d = {n, y, mue, beta};
    f.f = &curvefit_f;              // function f
    f.df = &curvefit_df;            // function df
    f.fdf = &curvefit_fdf;          // function fdf
    f.n = n;                        // observations
    f.p = p;                        // parameters
    f.params = &d;                  // struct
    
    T = gsl_multifit_fdfsolver_lmsder;                  // 
Levenberg-Marquardt-Solver type
    s = gsl_multifit_fdfsolver_alloc(T, n, p);          // 
Levenberg-Marquardt-Solver object
    gsl_vector_view x = gsl_vector_view_array(x_init, p);
    gsl_multifit_fdfsolver_set(s, &f, &x.vector);       // solver-collection
    curveFitPrintState(iter, s);                        // solving ... first run
    do {                                                // more runs 
        iter++;
        status = gsl_multifit_fdfsolver_iterate(s);
        printf("solver-status = %s\n", gsl_strerror(status));
        curveFitPrintState(iter, s);
        if (status) {
            break;
        } // break if minimum is achieved
        status = gsl_multifit_test_delta(s->dx, s->x, 1e-4, 1e-4);
    } while (status == GSL_CONTINUE && iter < 500);
    
    // Print Result of solver
    
    // Analyse Solving Results
    
    // free memory 
    gsl_multifit_fdfsolver_free(s);
    gsl_matrix_free(covar);
    gsl_vector_free( uvalsReal );
    gsl_rng_free(r);
    
    cout << " ... GSL-Curve-Fitting-End \n";
    return 0;
}

--
Jochen Bauer