Levenberg-Marquardt-Algorithmus - welche Implementierung für C++



  • Hi,

    eine Implementierung ist Teil der Gnu Scientific Library, aber man findet kaum etwas dazu. Außerdem kommen viele Funktionen vor, die an C erinnern, printf() und Speicherfreigaben. Ebenso wenn ich nach den spezifischen Funktionen der Bibliothek google, finde ich sehr wenig. Ich glaube also, es gibt eine andere verbreitetere Umsetzung für C++-Programme, bei denen ggf. mehr Tutorials und Doku zu finden ist. Danke für den ein oder anderen Tipp. Hier noch der Link zur GSL-Variante:
    http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Nonlinear-Least_002dSquares-Fitting.html


  • Mod

    Ich fürchte in der Numerik kommst du um C (oder gar Fortran!) Schnittstellen nicht so leicht herum. Wenn jemand viel Arbeit in eine allgemein anwendbare Bibliothek steckt, dann macht er das Interface so, dass es von möglichst vielen Sprachen benutzt werden kann und C ist meistens der kleinste gemeinsame Nenner. Ebenso natürlich die Beispiele.

    Meine Empfehlung wäre, das einfach so zu schlucken und sich des C-Interfaces zu bedienen (es hindert dich ja auch nichts, z.B. einen vector anstatt eines dynamischen Arrays zu nehmen, dann übergibst du eben die Adresse des ersten Elements als Parameter). Nur wenn du es wirklich an sehr vielen Stellen sehr oft brauchst, dann würde ich mir selber einen C++-Wrapper schreiben oder mal schauen ob es schon einen fertigen gibt. Mit diesem Stichwort findest du gewiss was bei Google.



  • Mein Tip:

    selber machen wenn es performant sein soll (Ist auch nichts anderes als Gauss-Newton). Sonst kannst du dir ja mal die C++ Implementierung von LevMar anschauen
    http://www.ics.forth.gr/~lourakis/levmar/



  • Shit, ich suche nun schon seit Tagen ... einmal fand ich etwas in Numerical Recipes, da ist es so, dass dort wiederum die Verteilung nicht implementiert ist. Außerdem ist die Variante langsam und es werden wohl sehr viele Berechnungen durchgeführt, es muss also schon schnell sein. In der Gnu Scientific Library-Variante ist sehr viel anzugeben und die Doku ist sehr gering, da blick ich kaum durch. Tutorials oder Beispiele zur GSL sind rar. Und levmar bekomme ich hier auf meiner Maschine (Ubuntu 10.04) nicht ans Laufen ... dazu muss der Code auf einer Maschine laufen, die ich nicht beeinflussen kann, also da kann ich nicht mal rumprobieren, um lange etwas zu installieren - geht mir das auf den Sack mittlerweile! Das gibts doch nicht, dass das hier noch keiner lösen musste oder dass das so schwierig ist.

    So - nächster Versuch: http://devernay.free.fr/hacks/cminpack/cminpack.html



  • Sodala,

    ich versuche mich gerade dran das GSL-Beispiel umzubiegen. Im Endeffekt kann das ja nicht so schwierig sein und da ich die GSL ohnehin einsetze. Ich habe die Daten im Code dazu einige Fragen eingestellt und die Links zu den Formeln der Gumbel-Verteilung mit den Parametern mue und beta. Vielleicht findet ja einer den Thread der vorher schon etwas ähnliches gemacht hat und kann mir den ein oder anderen Tipp geben. gslCurveFit() ist die main-Funktion.

    Hier ist der Link zum Beispiel und darunter mein Code:
    http://www.gnu.org/software/gsl/manual/html_node/Example-programs-for-Nonlinear-Least_002dSquares-Fitting.html

    // user-provided-functions: 
    // Gumbel-densitiy:     y = exp( ( -1 * ( x - mue ) ) / beta ) * exp( -1 * ( exp( ( -1 * ( x - mue ) ) / beta )) ) / beta ;
    // Gumbel-cumulative:   y = exp( - exp( - ( x - mue ) / beta ) ); 
    // Gumbel-Vals uniform-dice: y = mue - beta * ln( -1 * ln ( x ) );
    
    struct curvefitdata // data-container 
    {
        size_t n;           // observations
        double * y;         // y-val f(x)
        double * mue;       // mue, location of distribution
        double * beta;      // beta, scale of distribution 
    };
    
    // function stores f(x)-vals to incoming x-vals. 
    int curvefit_f( const gsl_vector * x, void *data,  gsl_vector * f ) 
    {
        // Attention: value x comes from parameter 'data', GSL-x represents variable to look for, often the y in math!
        // params: 
        //  x = current position, that means observation // TODO is this right?!
        //  f = current value at current position // lots of 
        //  data = current-data-container, see struct curvefitdata
    
        size_t n = ((struct curvefitdata *) data)->n;
        double *y = ((struct curvefitdata *) data)->y;
        double *mue = ((struct curvefitdata *) data)->mue;
        double *beta = ((struct curvefitdata *) data)->beta;
    
        // pick datacontainer x, the vector
        // seems this are the weights and above are the params to fit
        double w1 = gsl_vector_get(x, 0);
        double w2 = gsl_vector_get(x, 1);
    
        size_t i; // TODO Are those number of obervations?
        for (i = 0; i < n; i++) 
        {
            // Modell: 
            // y = exp( ( -1 * ( x - mue ) ) / beta ) * exp( -1 * ( exp( ( -1 * ( x - mue ) ) / beta )) ) / beta ;
    
            // TODO Check the meaning of this part of the function
            double t = i; // TODO is t for time and should this be implicit type casting?
            // double Yi = exp( ( -1 * ( x - mue ) ) / beta ) * exp( -1 * ( exp( ( -1 * ( x - mue ) ) / beta )) ) / beta ;
            double Yi = 3; // TODO change to formula
            gsl_vector_set(f, i, (Yi - y[i]) / beta[i]); // 
        }
        return GSL_SUCCESS;
    }
    
    // function stores n-by-p-values of matrix
    int curvefit_df( const gsl_vector * x, void *data,  gsl_matrix * J) 
    {   
        size_t n = ((struct curvefitdata *) data)->n;
        double *y = ((struct curvefitdata *) data)->y;
        double *mue = ((struct curvefitdata *) data)->mue;
        double *beta = ((struct curvefitdata *) data)->beta;
    
        // TODO see ref
        double A = gsl_vector_get(x, 0);
        double lambda = gsl_vector_get(x, 1);
        size_t i;
    
        for (i = 0; i < n; i++) {
    
            // TODO, see ref
    
            //double t = i;
            //double s = mue[i];
            //double e = exp(-lambda * t);
            //gsl_matrix_set(J, i, 0, e / s);
            //gsl_matrix_set(J, i, 1, -t * A * e / s);
            //gsl_matrix_set(J, i, 2, 1 / s);
        }
        return GSL_SUCCESS;
    }
    
    // function calls curvefit_f and curvefit_df simultanously
    int curvefit_fdf( const gsl_vector * x, void *data, gsl_vector * f, gsl_matrix * J) {
    
        curvefit_f(x, data, f);
        curvefit_df(x, data, J);
        return GSL_SUCCESS;
    }
    
    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));
    }
    
    void gslCurveFit()
    {
        // IDE-Properties:
        // Netbeans-Linker-Additional-Options: -lgsl -lgslcblas -lm
    
        // References:
        // 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
        cout << "GSL Curve Fitting-Example \n" << endl;
    
        // 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 ) );
        }
    
        // functions-components
        // struct curvefitdata - parameters and number of observations
        // 3 corefunctions x_f, x_df and x_fdf, got hooked to solver
        // solver attributes are
        //          vektor x current position, 
        //          vektor f current value at current position; 
        //          vektor dx difference after one step down to local minimum
        //          matrix J Jacobi-Matrix at current position
    
        // main() 
        // solver get initialized, data and starting guesses get to it as well.
    
        // Curvefitting to Gumbel-Dist, Levenberg-Marquardt
        int status;         // solver-status
        int i;              // index
        int iter = 0;       // counter for steps down to minimum
        int constn = 30;    // number of observations, 'N' in Documentation
        const size_t n = constn;    
        const size_t p = 2;         // parameters, mue and beta
        const gsl_multifit_fdfsolver_type * T = gsl_multifit_fdfsolver_lmder;       // solver-creation part 1
        gsl_multifit_fdfsolver * s = gsl_multifit_fdfsolver_alloc(T, constn, 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 functions
    
        // 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); // Solver verwaltet die Abweichungen
        } while (status == GSL_CONTINUE && iter < 500);
    
        // Print Result of solver
        gsl_multifit_covar (s->J, 0.0, covar);
    
            #define FIT(i) gsl_vector_get(s->x, i)
            #define ERR(i) sqrt(gsl_matrix_get(covar,i,i))
    
        {
            double chi = gsl_blas_dnrm2(s->f);
            double dof = n - p;
            double c = GSL_MAX_DBL(1, chi / sqrt(dof));
            printf("chisq/dof = %g\n", pow(chi, 2.0) / dof);
            printf("mue = % .5f + / - % .5f\n",     FIT(0),         c * ERR(0));
            printf("beta = %.5f +/- %.5f\n",        FIT(1),         c * ERR(1));
        }
        printf ("status = %s\n", gsl_strerror (status));
    
        // 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";
    }
    


  • Nächste Levenberg-Marquardt-Umsetzung, diese erfordert wirklich nur die Angabe einer Funktion die dann minimiert wird, mal schauen, ob ich damit klar komme: http://www.quantcode.com/modules/mydownloads/singlefile.php?lid=436



  • So habe es nun zusammengehämmert, da der Code sonst super-umfassend ist, habe ich Quelltextkommentare rausgenommen, bitte bedenken, einiges an Code ist von Fremprojekten. Es werden 1000 Werte aus einem Experiment eingelesen, daraus wird mittels GSL dann arithmetisches Mittel und Standardabweichung ermittelt. Es wird angenommen, dass es sich um eine Gumbel-Verteilung handelt. Daher können direkt mit den bekannten Größen die Parameter mue und beta berechnet werden. Okay, dann werden mittels GSL Histogramme angelegt und die Bins und Höhen kann man nutzen um die Parameter der kumulierten Verteilung via Quantlib-C-Minpack-Implementierung zu berechnen. Dazu kommt noch ein Test, ob die Gumbel-Verteilung zutrifft oder nicht, hier nehme ich den KS-Test aus Numerical Recipes, da gibts aber noch Problemchen. Da da heut ein Riesengefrickel war, poste ich hier mal den Code, vielleicht hilft die Vorgehensweise ja jemand.

    #include <algorithm>
    #include <cmath>
    #include <cstdio>
    #include <fstream>
    #include <gsl/gsl_blas.h>
    #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_sort.h>
    #include <gsl/gsl_statistics.h>
    #include <gsl/gsl_vector.h>
    #include <levmar.h>
    #include <math.h>
    #include <iostream>
    #include <sstream>
    #include <stdlib.h>  
    #include <stdio.h>
    #include <string.h>
    #include <time.h>
    #include <vector>
    #include "nr3.h"
    #include "sort.h"
    
    using namespace std;
    
    // Quantlib-globale Variablen
    double* independent_variables;
    double* oberved_values;
    
    // XXXXXXXXXXXXXXX Start Quantlib-Header lmdif XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    /************************lmdif*************************/
    namespace QuantLib {
      namespace MINPACK {
    #define BUG 0
    /* resolution of arithmetic */
    double MACHEP = 1.2e-16;
    /* smallest nonzero number */
    double DWARF = 1.0e-38;
    double enorm(int n,double* x) 
    {
                int i;
                double agiant, floatn, s1, s2, s3, xabs, x1max, x3max;
                double ans, temp;
                static double rdwarf = 3.834e-20;
                static double rgiant = 1.304e19;
                static double zero = 0.0;
                static double one = 1.0;
    s1 = zero;
    s2 = zero;
    s3 = zero;
    x1max = zero;
    x3max = zero;
    floatn = n;
    agiant = rgiant/floatn;
    for( i=0; i<n; i++ )
    {
    xabs = std::fabs(x[i]);
    if( (xabs > rdwarf) && (xabs < agiant) )
        {
        s2 += xabs*xabs;
        continue;
        }
    
    if(xabs > rdwarf)
        {
        if(xabs > x1max)
            {
            temp = x1max/xabs;
            s1 = one + s1*temp*temp;
            x1max = xabs;
            }
        else
            {
            temp = xabs/x1max;
            s1 += temp*temp;
            }
        continue;
        }
    
    if(xabs > x3max)
        {
        temp = x3max/xabs;
        s3 = one + s3*temp*temp;
        x3max = xabs;
        }
    else
        {
        if(xabs != zero)
            {
            temp = xabs/x3max;
            s3 += temp*temp;
            }
        }
    }
    if(s1 != zero)
        {
        temp = s1 + (s2/x1max)/x1max;
        ans = x1max*std::sqrt(temp);
        return(ans);
        }
    if(s2 != zero)
        {
        if(s2 >= x3max)
            temp = s2*(one+(x3max/s2)*(x3max*s3));
        else
            temp = x3max*((s2/x3max)+(x3max*s3));
        ans = std::sqrt(temp);
        }
    else
        {
        ans = x3max*std::sqrt(s3);
        }
    return(ans);
    }
    /************************lmmisc.c*************************/
    double dmax1(double a,double b)
    {
    if( a >= b )
        return(a);
    else
        return(b);
    }
    
    double dmin1(double a,double b)
    {
    if( a <= b )
        return(a);
    else
        return(b);
    }
    
    int min0(int a,int b)
    {
    if( a <= b )
        return(a);
    else
        return(b);
    }
    
    int mod( int k, int m )
    {
    return( k % m );
    }
    
    #if BUG
    void pmat( int m, int n, double* y  )
    {
    int i, j, k;
    
    k = 0;
    for( i=0; i<m; i++ )
        {
        for( j=0; j<n; j++ )
            {
            printf( "%.5e ", y[k] );
            k += 1;
            }
        printf( "\n" );
        }
    }
    #endif
    
    // XXXXXXX Fitting-Model XXXXXXXXXX
    /***********Sample of user supplied function****************
     * m = number of functions
     * n = number of variables
     * x = vector of function arguments
     * fvec = vector of function values
     * iflag = error return variable
     */
    void fcn(int m,int n, double* x, double* fvec,int *iflag)
    {
        double a=x[0]; 
        double b=x[1]; 
        for (int i=0;i<m;i++)
        {
            double xi=independent_variables[i];
            double observed_value=oberved_values[i];
            //double fitted_value=a*xi*xi*xi + b*xi*xi + c; // Tut-Code
            //double ydist = exp( - exp( - ( x - gumbelmue ) / gumbelbeta ) ); // cumulated
            //double ydist = exp( ( -1 * ( x - gumbelmue ) ) / gumbelbeta ) * exp( -1 * ( exp( ( -1 * ( x - gumbelmue ) ) / gumbelbeta )) ) / gumbelbeta ;
            //double gumbelDouble = ydist * ( elements / bins );
            double fitted_value=exp( - exp( - ( xi - a ) / b ) );
            fvec[i]=fitted_value-observed_value;
        }
    
    }
    void
    fdjac2(int m,int n,double* x,double* fvec,double* fjac,int,
           int* iflag,double epsfcn,double* wa)
    {
    int i,j,ij;
    double eps,h,temp;
    static double zero = 0.0;
    temp = dmax1(epsfcn,MACHEP);
    eps = std::sqrt(temp);
    #if BUG
    printf( "fdjac2\n" );
    #endif
    ij = 0;
    for( j=0; j<n; j++ )
        {
        temp = x[j];
        h = eps * std::fabs(temp);
        if(h == zero)
            h = eps;
        x[j] = temp + h;
        fcn(m,n,x,wa,iflag);
        if( *iflag < 0)
            return;
        x[j] = temp;
        for( i=0; i<m; i++ )
            {
            fjac[ij] = (wa[i] - fvec[i])/h;
            ij += 1;    /* fjac[i+m*j] */
            }
        }
    #if BUG
    pmat( m, n, fjac );
    #endif
    }
    /************************qrfac.c*************************/
    void
    qrfac(int m,int n,double* a,int,int pivot,int* ipvt,
          int,double* rdiag,double* acnorm,double* wa)
    {
    int i,ij,jj,j,jp1,k,kmax,minmn;
    double ajnorm,sum,temp;
    static double zero = 0.0;
    static double one = 1.0;
    static double p05 = 0.05;
    ij = 0;
    for( j=0; j<n; j++ )
        {
        acnorm[j] = enorm(m,&a[ij]);
        rdiag[j] = acnorm[j];
        wa[j] = rdiag[j];
        if(pivot != 0)
            ipvt[j] = j;
        ij += m; /* m*j */
        }
    #if BUG
    printf( "qrfac\n" );
    #endif
    
    minmn = min0(m,n);
    for( j=0; j<minmn; j++ )
    {
    if(pivot == 0)
        goto L40;
    
    kmax = j;
    for( k=j; k<n; k++ )
        {
        if(rdiag[k] > rdiag[kmax])
            kmax = k;
        }
    if(kmax == j)
        goto L40;
    
    ij = m * j;
    jj = m * kmax;
    for( i=0; i<m; i++ )
        {
        temp = a[ij]; /* [i+m*j] */
        a[ij] = a[jj]; /* [i+m*kmax] */
        a[jj] = temp;
        ij += 1;
        jj += 1;
        }
    rdiag[kmax] = rdiag[j];
    wa[kmax] = wa[j];
    k = ipvt[j];
    ipvt[j] = ipvt[kmax];
    ipvt[kmax] = k;
    L40:
    jj = j + m*j;
    ajnorm = enorm(m-j,&a[jj]);
    if(ajnorm == zero)
        goto L100;
    if(a[jj] < zero)
        ajnorm = -ajnorm;
    ij = jj;
    for( i=j; i<m; i++ )
        {
        a[ij] /= ajnorm;
        ij += 1; /* [i+m*j] */
        }
    a[jj] += one;
    
    jp1 = j + 1;
    if(jp1 < n )
    {
    for( k=jp1; k<n; k++ )
        {
        sum = zero;
        ij = j + m*k;
        jj = j + m*j;
        for( i=j; i<m; i++ )
            {
            sum += a[jj]*a[ij];
            ij += 1; /* [i+m*k] */
            jj += 1; /* [i+m*j] */
            }
        temp = sum/a[j+m*j];
        ij = j + m*k;
        jj = j + m*j;
        for( i=j; i<m; i++ )
            {
            a[ij] -= temp*a[jj];
            ij += 1; /* [i+m*k] */
            jj += 1; /* [i+m*j] */
            }
        if( (pivot != 0) && (rdiag[k] != zero) )
            {
            temp = a[j+m*k]/rdiag[k];
            temp = dmax1( zero, one-temp*temp );
            rdiag[k] *= std::sqrt(temp);
            temp = rdiag[k]/wa[k];
            if( (p05*temp*temp) <= MACHEP)
                {
                rdiag[k] = enorm(m-j-1,&a[jp1+m*k]);
                wa[k] = rdiag[k];
                }
            }
        }
    }
    L100:
        rdiag[j] = -ajnorm;
    }
    }
    /************************qrsolv.c*************************/
    void
    qrsolv(int n,double* r,int ldr,int* ipvt,double* diag,double* qtb,
           double* x,double* sdiag,double* wa)
    {
    int i,ij,ik,kk,j,jp1,k,kp1,l,nsing;
    double cos,cotan,qtbpj,sin,sum,tan,temp;
    static double zero = 0.0;
    static double p25 = 0.25;
    static double p5 = 0.5;
    kk = 0;
    for( j=0; j<n; j++ )
        {
        ij = kk;
        ik = kk;
        for( i=j; i<n; i++ )
            {
            r[ij] = r[ik];
            ij += 1;   /* [i+ldr*j] */
            ik += ldr; /* [j+ldr*i] */
            }
        x[j] = r[kk];
        wa[j] = qtb[j];
        kk += ldr+1; /* j+ldr*j */
        }
    #if BUG
    printf( "qrsolv\n" );
    #endif
    for( j=0; j<n; j++ )
    {
    l = ipvt[j];
    if(diag[l] == zero)
        goto L90;
    for( k=j; k<n; k++ )
        sdiag[k] = zero;
    sdiag[j] = diag[l];
    qtbpj = zero;
    for( k=j; k<n; k++ )
        {
        if(sdiag[k] == zero)
            continue;
        kk = k + ldr * k;
        if(std::fabs(r[kk]) < std::fabs(sdiag[k]))
            {
            cotan = r[kk]/sdiag[k];
            sin = p5/std::sqrt(p25+p25*cotan*cotan);
            cos = sin*cotan;
            }
        else
            {
            tan = sdiag[k]/r[kk];
            cos = p5/std::sqrt(p25+p25*tan*tan);
            sin = cos*tan;
            }
        r[kk] = cos*r[kk] + sin*sdiag[k];
        temp = cos*wa[k] + sin*qtbpj;
        qtbpj = -sin*wa[k] + cos*qtbpj;
        wa[k] = temp;
        kp1 = k + 1;
        if( n > kp1 )
            {
            ik = kk + 1;
            for( i=kp1; i<n; i++ )
                {
                temp = cos*r[ik] + sin*sdiag[i];
                sdiag[i] = -sin*r[ik] + cos*sdiag[i];
                r[ik] = temp;
                ik += 1; /* [i+ldr*k] */
                }
            }
        }
    L90:
        kk = j + ldr*j;
        sdiag[j] = r[kk];
        r[kk] = x[j];
    }
    nsing = n;
    for( j=0; j<n; j++ )
        {
        if( (sdiag[j] == zero) && (nsing == n) )
            nsing = j;
        if(nsing < n)
            wa[j] = zero;
        }
    if(nsing < 1)
        goto L150;
    
    for( k=0; k<nsing; k++ )
        {
        j = nsing - k - 1;
        sum = zero;
        jp1 = j + 1;
        if(nsing > jp1)
            {
            ij = jp1 + ldr * j;
            for( i=jp1; i<nsing; i++ )
                {
                sum += r[ij]*wa[i];
                ij += 1; /* [i+ldr*j] */
                }
            }
        wa[j] = (wa[j] - sum)/sdiag[j];
        }
    L150:
    for( j=0; j<n; j++ )
        {
        l = ipvt[j];
        x[l] = wa[j];
        }
    }
    /************************lmpar.c*************************/
    void
    lmpar(int n,double* r,int ldr,int* ipvt,double* diag,
          double* qtb,double delta,double* par,double* x,double* sdiag,
          double* wa1,double* wa2)
    {
    int i,iter,ij,jj,j,jm1,jp1,k,l,nsing;
    double dxnorm,fp,gnorm,parc,parl,paru;
    double sum,temp;
    static double zero = 0.0;
    static double p1 = 0.1;
    static double p001 = 0.001;
    extern double DWARF;
    #if BUG
    printf( "lmpar\n" );
    #endif
    nsing = n;
    jj = 0;
    for( j=0; j<n; j++ )
        {
        wa1[j] = qtb[j];
        if( (r[jj] == zero) && (nsing == n) )
            nsing = j;
        if(nsing < n)
            wa1[j] = zero;
        jj += ldr+1; /* [j+ldr*j] */
        }
    #if BUG
    printf( "nsing %d ", nsing );
    #endif
    if(nsing >= 1)
        {
        for( k=0; k<nsing; k++ )
            {
            j = nsing - k - 1;
            wa1[j] = wa1[j]/r[j+ldr*j];
            temp = wa1[j];
            jm1 = j - 1;
            if(jm1 >= 0)
                {
                ij = ldr * j;
                for( i=0; i<=jm1; i++ )
                    {
                    wa1[i] -= r[ij]*temp;
                    ij += 1;
                    }
                }
            }
        }
    for( j=0; j<n; j++ )
        {
        l = ipvt[j];
        x[l] = wa1[j];
        }
    iter = 0;
    for( j=0; j<n; j++ )
        wa2[j] = diag[j]*x[j];
    dxnorm = enorm(n,wa2);
    fp = dxnorm - delta;
    if(fp <= p1*delta)
        {
    #if BUG
        printf( "going to L220\n" );
    #endif
        goto L220;
        }
    parl = zero;
    if(nsing >= n)
        {
        for( j=0; j<n; j++ )
            {
            l = ipvt[j];
            wa1[j] = diag[l]*(wa2[l]/dxnorm);
            }
        jj = 0;
        for( j=0; j<n; j++ )
            {
            sum = zero;
            jm1 = j - 1;
            if(jm1 >= 0)
                {
                ij = jj;
                for( i=0; i<=jm1; i++ )
                    {
                    sum += r[ij]*wa1[i];
                    ij += 1;
                    }
                }
            wa1[j] = (wa1[j] - sum)/r[j+ldr*j];
            jj += ldr; /* [i+ldr*j] */
            }
        temp = enorm(n,wa1);
        parl = ((fp/delta)/temp)/temp;
        }
    jj = 0;
    for( j=0; j<n; j++ )
        {
        sum = zero;
        ij = jj;
        for( i=0; i<=j; i++ )
            {
            sum += r[ij]*qtb[i];
            ij += 1;
            }
        l = ipvt[j];
        wa1[j] = sum/diag[l];
        jj += ldr; /* [i+ldr*j] */
        }
    gnorm = enorm(n,wa1);
    paru = gnorm/delta;
    if(paru == zero)
        paru = DWARF/dmin1(delta,p1);
    *par = dmax1( *par,parl);
    *par = dmin1( *par,paru);
    if( *par == zero)
        *par = gnorm/dxnorm;
    #if BUG
    printf( "parl %.4e  par %.4e  paru %.4e\n", parl, *par, paru );
    #endif
    L150:
    iter += 1;
    if( *par == zero)
        *par = dmax1(DWARF,p001*paru);
    temp = std::sqrt( *par );
    for( j=0; j<n; j++ )
        wa1[j] = temp*diag[j];
    qrsolv(n,r,ldr,ipvt,wa1,qtb,x,sdiag,wa2);
    for( j=0; j<n; j++ )
        wa2[j] = diag[j]*x[j];
    dxnorm = enorm(n,wa2);
    temp = fp;
    fp = dxnorm - delta;
    
    if( (std::fabs(fp) <= p1*delta)
     || ((parl == zero) && (fp <= temp) && (temp < zero))
     || (iter == 10) )
        goto L220;
    
    for( j=0; j<n; j++ )
        {
        l = ipvt[j];
        wa1[j] = diag[l]*(wa2[l]/dxnorm);
        }
    jj = 0;
    for( j=0; j<n; j++ )
        {
        wa1[j] = wa1[j]/sdiag[j];
        temp = wa1[j];
        jp1 = j + 1;
        if(jp1 < n)
            {
            ij = jp1 + jj;
            for( i=jp1; i<n; i++ )
                {
                wa1[i] -= r[ij]*temp;
                ij += 1; /* [i+ldr*j] */
                }
            }
        jj += ldr; /* ldr*j */
        }
    temp = enorm(n,wa1);
    parc = ((fp/delta)/temp)/temp;
    if(fp > zero)
        parl = dmax1(parl, *par);
    if(fp < zero)
        paru = dmin1(paru, *par);
    *par = dmax1(parl, *par + parc);
    goto L150;
    L220:
    if(iter == 0)
        *par = zero;
    }
    /*********************** lmdif.c ****************************/
    
    void lmdif(int m,int n,double* x,double* fvec,double ftol,
          double xtol,double gtol,int maxfev,double epsfcn,
          double* diag, int mode, double factor,
          int nprint, int* info,int* nfev,double* fjac,
          int ldfjac,int* ipvt,double* qtf,
          double* wa1,double* wa2,double* wa3,double* wa4)
    {
    int i,iflag,ij,jj,iter,j,l;
    double actred,delta=0,dirder,fnorm,fnorm1,gnorm;
    double par,pnorm,prered,ratio;
    double sum,temp,temp1,temp2,temp3,xnorm=0;
    static double one = 1.0;
    static double p1 = 0.1;
    static double p5 = 0.5;
    static double p25 = 0.25;
    static double p75 = 0.75;
    static double p0001 = 1.0e-4;
    static double zero = 0.0;
    *info = 0;
    iflag = 0;
    *nfev = 0;
    
    if( (n <= 0) || (m < n) || (ldfjac < m) || (ftol < zero)
        || (xtol < zero) || (gtol < zero) || (maxfev <= 0)
        || (factor <= zero) )
        goto L300;
    
    if( mode == 2 )
        { /* scaling by diag[] */
        for( j=0; j<n; j++ )
            {
            if( diag[j] <= 0.0 )
                goto L300;
            }
        }
    #if BUG
    printf( "lmdif\n" );
    #endif
    iflag = 1;
    fcn(m,n,x,fvec,&iflag);
    *nfev = 1;
    if(iflag < 0)
        goto L300;
    fnorm = enorm(m,fvec);
    par = zero;
    iter = 1;
    L30:
    iflag = 2;
    fdjac2(m,n,x,fvec,fjac,ldfjac,&iflag,epsfcn,wa4);
    *nfev += n;
    if(iflag < 0)
        goto L300;
    if( nprint > 0 )
        {
        iflag = 0;
        if(mod(iter-1,nprint) == 0)
            {
            fcn(m,n,x,fvec,&iflag);
            if(iflag < 0)
                goto L300;
            #if BUG
            printf( "fnorm %.15e\n", enorm(m,fvec) );
            #endif
            }
        }
    qrfac(m,n,fjac,ldfjac,1,ipvt,n,wa1,wa2,wa3);
    
    if(iter == 1)
        {
        if(mode != 2)
            {
            for( j=0; j<n; j++ )
                {
                diag[j] = wa2[j];
                if( wa2[j] == zero )
                    diag[j] = one;
                }
            }
        for( j=0; j<n; j++ )
            wa3[j] = diag[j] * x[j];
    
        xnorm = enorm(n,wa3);
        delta = factor*xnorm;
        if(delta == zero)
            delta = factor;
        }
    
    for( i=0; i<m; i++ )
        wa4[i] = fvec[i];
    jj = 0;
    for( j=0; j<n; j++ )
        {
        temp3 = fjac[jj];
        if(temp3 != zero)
            {
            sum = zero;
            ij = jj;
            for( i=j; i<m; i++ )
                {
                sum += fjac[ij] * wa4[i];
                ij += 1;    /* fjac[i+m*j] */
                }
            temp = -sum / temp3;
            ij = jj;
            for( i=j; i<m; i++ )
                {
                wa4[i] += fjac[ij] * temp;
                ij += 1;    /* fjac[i+m*j] */
                }
            }
        fjac[jj] = wa1[j];
        jj += m+1;  /* fjac[j+m*j] */
        qtf[j] = wa4[j];
        }
     gnorm = zero;
     if(fnorm != zero)
        {
        jj = 0;
        for( j=0; j<n; j++ )
            {
            l = ipvt[j];
            if(wa2[l] != zero)
                {
                sum = zero;
                ij = jj;
                for( i=0; i<=j; i++ )
                    {
                    sum += fjac[ij]*(qtf[i]/fnorm);
                    ij += 1; /* fjac[i+m*j] */
                    }
                gnorm = dmax1(gnorm,std::fabs(sum/wa2[l]));
                }
            jj += m;
            }
        }
     if(gnorm <= gtol)
        *info = 4;
     if( *info != 0)
        goto L300;
    
     if(mode != 2)
        {
        for( j=0; j<n; j++ )
            diag[j] = dmax1(diag[j],wa2[j]);
        }
    L200:
    lmpar(n,fjac,ldfjac,ipvt,diag,qtf,delta,&par,wa1,wa2,wa3,wa4);
    
    for( j=0; j<n; j++ )
        {
           wa1[j] = -wa1[j];
           wa2[j] = x[j] + wa1[j];
           wa3[j] = diag[j]*wa1[j];
        }
    pnorm = enorm(n,wa3);
    
    if(iter == 1)
        delta = dmin1(delta,pnorm);
    
    iflag = 1;
    fcn(m,n,wa2,wa4,&iflag);
    *nfev += 1;
    if(iflag < 0)
        goto L300;
    fnorm1 = enorm(m,wa4);
    #if BUG
    printf( "pnorm %.10e  fnorm1 %.10e\n", pnorm, fnorm1 );
    #endif
    
    actred = -one;
    if( (p1*fnorm1) < fnorm)
        {
        temp = fnorm1/fnorm;
        actred = one - temp * temp;
        }
    
    jj = 0;
    for( j=0; j<n; j++ )
        {
        wa3[j] = zero;
        l = ipvt[j];
        temp = wa1[l];
        ij = jj;
        for( i=0; i<=j; i++ )
            {
            wa3[i] += fjac[ij]*temp;
            ij += 1; /* fjac[i+m*j] */
            }
        jj += m;
        }
    temp1 = enorm(n,wa3)/fnorm;
    temp2 = (std::sqrt(par)*pnorm)/fnorm;
    prered = temp1*temp1 + (temp2*temp2)/p5;
    dirder = -(temp1*temp1 + temp2*temp2);
    
    ratio = zero;
    if(prered != zero)
        ratio = actred/prered;
    
    if(ratio <= p25)
        {
        if(actred >= zero)
            temp = p5;
        else
            temp = p5*dirder/(dirder + p5*actred);
        if( ((p1*fnorm1) >= fnorm)
        || (temp < p1) )
            temp = p1;
           delta = temp*dmin1(delta,pnorm/p1);
           par = par/temp;
        }
    else
        {
        if( (par == zero) || (ratio >= p75) )
            {
            delta = pnorm/p5;
            par = p5*par;
            }
        }
    
    if(ratio >= p0001)
        {
        for( j=0; j<n; j++ )
            {
            x[j] = wa2[j];
            wa2[j] = diag[j]*x[j];
            }
        for( i=0; i<m; i++ )
            fvec[i] = wa4[i];
        xnorm = enorm(n,wa2);
        fnorm = fnorm1;
        iter += 1;
        }
    if( (std::fabs(actred) <= ftol)
      && (prered <= ftol)
      && (p5*ratio <= one) )
        *info = 1;
    if(delta <= xtol*xnorm)
        *info = 2;
    if( (std::fabs(actred) <= ftol)
      && (prered <= ftol)
      && (p5*ratio <= one)
      && ( *info == 2) )
        *info = 3;
    if( *info != 0)
        goto L300;
    if( *nfev >= maxfev)
        *info = 5;
    if( (std::fabs(actred) <= MACHEP)
      && (prered <= MACHEP)
      && (p5*ratio <= one) )
        *info = 6;
    if(delta <= MACHEP*xnorm)
        *info = 7;
    if(gnorm <= MACHEP)
        *info = 8;
    if( *info != 0)
        goto L300;
    if(ratio < p0001)
        goto L200;
    goto L30;
    L300:
    if(iflag < 0)
        *info = iflag;
    iflag = 0;
    if(nprint > 0)
        fcn(m,n,x,fvec,&iflag);
    }}}
    // XXXXXXXXXXXXXXXXXXXXXXXXXX ENDE Quantlib-Header lmdif XXXXXXXXXXXXXXXXXXXXXXX
    
    // XXXXXXXXXXXXXXXXXXXXX Start Numerical Recipe KS-Test XXXXXXXXXXXXXXXXXXXXXXXX
    // NR-Struktur fuer KS-Verteilung
    struct KSdist {
        // Kap 14.3.3. S.737 + Kap 6.14.56 S. 334 
    
        double pks(double z) {
            // Return cumulative distribution function.
            if (z < 0.) throw ("bad z in KSdist");
            if (z == 0.) return 0.;
            if (z < 1.18) {
                double y = exp(-1.23370055013616983 / SQR(z));
                return 2.25675833419102515 * sqrt(-log(y))*(y + pow(y, 9) + pow(y, 25) + pow(y, 49));
            } else {
                double x = exp(-2. * SQR(z));
                return 1. - 2. * (x - pow(x, 4) + pow(x, 9));
            }
        }
    
        // Return complementary cumulative distribution function.
    
        double qks(double z) {
            if (z < 0.) throw ("bad z in KSdist");
            if (z == 0.) return 1.;
            if (z < 1.18) return 1. - pks(z);
            double x = exp(-2. * (z * z));
            return 2. * (x - pow(x, 4) + pow(x, 9));
        }
    
        // more code in the book ...
    };
    
    // NR-Funkton zum KS-Test
    void kstwo(VecDoub_IO &data1, VecDoub_IO &data2, Doub &d, Doub &prob) {
        // Given an array data1[0..n1-1], and an array data2[0..n2-1], this routine returns the K–S
        // statistic d and the p-value prob for the null hypothesis that the data sets are drawn from the
        // same distribution. Small values of prob show that the cumulative distribution function of data1
        // is significantly different from that of data2. The arrays data1 and data2 are modified by being
        // sorted into ascending order.
        int j1 = 0, j2 = 0, n1 = data1.size(), n2 = data2.size();
        double d1, d2, dt, en1, en2, en, fn1 = 0.0, fn2 = 0.0;
        KSdist ks;
        sort(data1);
        sort(data2);
        // Ausgabe der beiden sortierten Collections
        //int nrRealVecLength = data1.size();
        //int nrGumbelVecLength = data2.size();
        //for (int i = 0; i < nrRealVecLength; i++) {
        //   cout << "Debug nrRealVec data1 ... Lauf " << i << ": " << data1[i] << endl;
        //}
        //for (int i = 0; i < nrGumbelVecLength; i++) {
        //    cout << "Debug nrGumbelVec data2 ... Lauf " << i << ": " << data2[i] << endl;
        //}
    
        en1 = n1;
        en2 = n2;
        d = 0.0;
        while (j1 < n1 && j2 < n2) // If we are not done...
        {
            if ((d1 = data1[j1]) <= (d2 = data2[j2])) // Next step is in data1.
            {
                do {
                    fn1 = ++j1 / en1;
                } while (j1 < n1 && d1 == data1[j1]);
            }
            if (d2 <= d1) // Next step is in data2.
            {
                do {
                    fn2 = ++j2 / en2;
                } while (j2 < n2 && d2 == data2[j2]);
            }
            if ((dt = abs(fn2 - fn1)) > d) {
                //cout << " fn2 " << fn2 << "; fn1 " << fn1 << "; d " << d << endl; 
                d = dt;
            }
        }
        en = sqrt(en1 * en2 / (en1 + en2));
        //cout << "* en " << en << " aus en1 " << en1 << " und en2 " << en2 << endl;
        prob = ks.qks( (en + 0.12 + 0.11 / en) * d); // Compute p-value.
    }
    // XXXXXXXXXXXXXXXXX KS-Test Ende XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX
    
    double doubleFromString( const std::string& s ) 
    {
        std::istringstream stream(s);
        double t;
        stream >> t;
        return t;
    }
    
    double roundDouble(double number, unsigned int digits) {
        number *= pow(10, digits);
        if (number >= 0)
            number = floor(number + 0.5);
        else
            number = ceil(number - 0.5); 
        number /= pow(10, digits);
        return number;
    }
    
    int fitAndTest()
    {
        // References:
        // Curve-Fitting: C++-Minpack of Quantlib: quantcode.com/modules/mydownloads/singlefile.php?lid=436
        // Gumbel-Distribution: en.wikipedia.org/wiki/Gumbel_distribution
        // Histograms: gnu.org/software/gsl/manual/html_node/index.html#Top
        // Kolmogorov-Smirnov-Test: Numerical Recipes - kstwo()
        cout << "Fit And Test \n" << endl;
        vector<double> uvalsRealVector;     // data collection
        double elements = 1000;             // experimental data
        double bins = 100;                  // histogram-bins    
        fstream fstrm; // read data from file
        string path = "/home/joba/data01.txt";
        fstrm.open(path.c_str(), ::std::ios_base::in | ::std::ios_base::binary);
        if (!fstrm) {
            cout << "XXX File not found! \n";
            exit(0);
        } 
        string currentValString;
        while (!fstrm.eof()) 
        {
            getline(fstrm, currentValString); //cout << "Zeile : " << curuval << endl;
            if ( currentValString == "" ){ }
            else // add value to collection
            {
                double currentVal = doubleFromString( currentValString );
                uvalsRealVector.push_back( currentVal );
            }
        }
        fstrm.close();
    
        // prepare data for GSL
        int numberOfElements = uvalsRealVector.size(); 
        double uvalsRealDataArray[ numberOfElements ];
        for( int i = 0; i < numberOfElements; i++ )
        {
            uvalsRealDataArray[i] = uvalsRealVector.at(i);
        }
    
        // GSL, compute mean and sdev, gsl_stats_mean(dataarray, stepwidth, numberOfElements);
        double uvalsRealMean   = gsl_stats_mean(uvalsRealDataArray, 1, numberOfElements); 
        double uvalsRealSdev   = gsl_stats_sd(uvalsRealDataArray, 1, numberOfElements);
    
        // calculate mue and beta from mean and sd - see gumbel distribution
        double beta = ( ( uvalsRealSdev * sqrt(6) ) / M_PI );
        double mue = ( uvalsRealMean - ( M_EULER * beta ) );
        cout << " ... mue: " << mue << "; beta: " << beta << " (calculated via mean and sd);" << endl;
    
        // start
        double gumbelmue    = mue;                  // 0.0939307 bei den genauen 30, aus R
        double gumbelbeta   = beta;                 // 0.0231414 bei den genauen 30, aus R
    
        // ----
        // histogram - real data
        int numberofbins = bins;
        gsl_histogram * h = gsl_histogram_alloc(numberofbins); // hist-pointer, n bins, n is rangeElements -1 
        double range[(numberofbins + 1)];
        double currentBinElement = 0.0;
        for (int i = 0; i < (numberofbins + 1); i++) 
        {
            range[i] = currentBinElement;
            currentBinElement = currentBinElement + 0.01; // Case: 100 bins TODO fix to relative
        } // bin-width, lower border inkl., upper border exkl. 
        gsl_histogram_set_ranges(h, range, (numberofbins + 1));
        for (int i = 0; i < numberOfElements; i++) // push vals to bins
        {
            double uvalRandom = uvalsRealDataArray[i];
            gsl_histogram_increment(h, uvalRandom);
        }
        FILE * pFile; // write hist to file 
        string fileName = "/home/joba/curhistoreal.txt";
        pFile = fopen(fileName.c_str(), "w");
        gsl_histogram_fprintf(pFile, h, "%g", "%g"); // paras: FILE, const gsl_histogram * h, const char * range_format, const char * bin_format
        fclose(pFile);
        //cout << " ... Histogram realData is written! \n"; 
        // Betrachten via awk '{print $1, $3 ; print $2, $3}' /home/joba/curhisto.txt | graph -T X 
        // or via cat /home/joba/curhisto.txt
    
        // ----
        // cumulated histogram from noncumulated data, real data
        int numberofbinscum = bins;
        gsl_histogram * hcum = gsl_histogram_alloc(numberofbinscum); 
        double rangecum[(numberofbinscum + 1)];
        double currentBinElementCum = 0.0;
        for (int i = 0; i < (numberofbinscum + 1); i++) {
            rangecum[i] = currentBinElementCum;
            currentBinElementCum = currentBinElementCum + 0.01; // Case: 100 bins TODO fix to relative
        }
        gsl_histogram_set_ranges(hcum, rangecum, (numberofbinscum + 1));
        for (int i = 0; i < numberOfElements; i++) 
        {
            double uvalRandom = uvalsRealDataArray[i];
            gsl_histogram_increment(hcum, uvalRandom); 
        }
        double curCountsRealCum = 0.0; // summing up bins
        for (int i = 0; i < numberofbinscum; i++)
        {
            double lower = 0.0; // get mid of bin
            double upper = 0.0;
            gsl_histogram_get_range( hcum, i, &lower, &upper );
            double classmid = (upper + lower) / 2; 
            double curBinHits = gsl_histogram_get( hcum, i ); 
            gsl_histogram_accumulate( hcum, classmid, curCountsRealCum );
            curCountsRealCum = curCountsRealCum + curBinHits;
            //double curHits = gsl_histogram_get( hcum, i ); // TODO Debug-Only
            //cout << "Curhistocum - Klassenmitte " << classmid << ", Hits: " << curHits << endl;
        }
        FILE * pFileCum; 
        string fileNameCum = "/home/joba/curhistocum.txt";
        pFileCum = fopen(fileNameCum.c_str(), "w");
        gsl_histogram_fprintf(pFileCum, hcum, "%g", "%g"); 
        fclose(pFileCum);
        //cout << " ... Histogram realcum is written! \n";
    
        // ---
        // grab bin mids and bin counts of one or both histograms
        // maybe compute theoretical gumbel-height to bin mids
        // -> base for Minpack-Curve-Fit and KS-Test is available 
        // afterwards.
        int elementsToCompare = 100;                // number of bins - elements to compare
        vector<double> binMidList;
        VecDoub nrRealVec( elementsToCompare );     // data points of experiment
        VecDoub nrGumbelVec( elementsToCompare );   // data points of theoretical distribution 
        VecDoub nrRealVecCumulated( elementsToCompare );     // data points of experiment
        //VecDoub nrGumbelVecCumulated( elementsToCompare );   // data points of theoretical distribution 
        for ( int i = 0; i < elementsToCompare; i++ )
        {
            // TODO there are useless computation for example redundant mid-calcs
            // grab bin-mids of real data for KS-Test
            double currentRealClassLower = 0.0;
            double currentRealClassUpper = 0.0;
            gsl_histogram_get_range( h, i, &currentRealClassLower, &currentRealClassUpper );
            double currentRealClassMid = (currentRealClassUpper + currentRealClassLower) / 2;
            binMidList.push_back( currentRealClassMid ); // for Minpack-Fit
            double absBinElements = gsl_histogram_get( h, i); // bin-count
            nrRealVec[i] = absBinElements; // / relBinsElements;
    
            // compute-y-val of gumbel-distribution for KS-test
            double x = currentRealClassMid;
            //double ydist = exp( - exp( - ( x - gumbelmue ) / gumbelbeta ) ); // cumulated gumbel dist
            double ydist = exp( ( -1 * ( x - gumbelmue ) ) / gumbelbeta ) * exp( -1 * ( exp( ( -1 * ( x - gumbelmue ) ) / gumbelbeta )) ) / gumbelbeta ; // gumbel-dist density
            double gumbelDouble = ydist * ( elements / bins );
            nrGumbelVec[i] = roundDouble(gumbelDouble, 0); 
            //cout << " -> Run - " << i << "; x - " << currentRealClassMid << "; real - " << nrRealVec[i] << "; gumbel - "  << nrGumbelVec[i] << "; y-func: "  << ydist << "; y-binelerel-func: " << endl;
    
            // grab bin-mids of real cumulated data for Minpack-Fit
            double currentRealClassLowerCumulated = 0.0;
            double currentRealClassUpperCumulated = 0.0;
            gsl_histogram_get_range( hcum, i, &currentRealClassLowerCumulated, &currentRealClassUpperCumulated );
            double currentRealClassMidCumulated = (currentRealClassUpperCumulated + currentRealClassLowerCumulated) / 2;
            double absBinElementsCumulated = gsl_histogram_get( hcum, i); // bin-count
            nrRealVecCumulated[i] = absBinElementsCumulated; // / relBinsElements;
    
            // compute-y-val of gumbel-cum-distribution for Minpack-Fit
            //double xCumulated = currentRealClassMidCumulated;
            //double ydistCumulated = exp( - exp( - ( xCumulated - gumbelmue ) / gumbelbeta ) ); // cumulated gumbel dist
            //double ydist = exp( ( -1 * ( x - gumbelmue ) ) / gumbelbeta ) * exp( -1 * ( exp( ( -1 * ( x - gumbelmue ) ) / gumbelbeta )) ) / gumbelbeta ; // gumbel-dist density
            //double gumbelDoubleCumulated = ydistCumulated * ( elements / bins );
            //nrGumbelVecCumulated[i] = roundDouble(gumbelDoubleCumulated, 5); 
            //cout << " -> Run - " << i << "; xCumulated - " << currentRealClassMidCumulated << "; realCumulated - " << nrRealVecCumulated[i] << ";" << endl; // gumbelCumulated - "  << nrGumbelVecCumulated[i] << "; y-func: "  << ydistCumulated << "; y-binelerel-func: " << gumbelDoubleCumulated << endl;
    
        }
    
        // free histograms
        gsl_histogram_free( h );
        gsl_histogram_free( hcum );
    
        // KS-test on gumbel with mue and beta from calculation
        double nrDis        = 0.0;
        double nrProb       = 0.0;
        //cout << "V1 realData ... " << endl;
        //for ( int i = 0; i < nrRealVec.size(); i++ )
        //{
        //    cout << nrRealVec[i] << ", ";
        //}
        //cout << endl;
        //cout << "V1 gumbelData ... " << endl;
        //for ( int i = 0; i < nrGumbelVec.size(); i++ )
        //{
        //    cout << nrGumbelVec[i] << ", "; 
        //}
        //cout << endl;
    
        kstwo( nrRealVec, nrGumbelVec, nrDis, nrProb );
        cout << "KS-Test: Dist: " << nrDis << "; P: " << nrProb << endl;
    
        // ---
        // fit with Quantlib-Minpack to Gumbel-Distribution, take cumulated relative bin heights and cumulated model
        // model goes to fcn() in header-section of this file
        // Data to fit
        int m = 100;        // number of observations or functions
        int n = 2;          // number of parameters to fit, here mue and beta
        independent_variables = new double[m];      // this is the vector of independent variables, the x-values, bin-mids
        oberved_values = new double[m];             // this is the vector of observed variables to be fitted, the y-values of observations
    
        for( int i = 0; i < m; i++ )
        { 
            // set x and y value, x bin mid, y = realCumBinHeight
            independent_variables[i] = binMidList.at( i );
            oberved_values[i] = nrRealVecCumulated[ i ] / 1000.0; // TODO fix to variable
        }
        double* x=new double[n]; // initial estimate of parameters vector
        x[0]=0.01;          // mue-Guess
        x[1]=0.005;         // beta-Guess
        double* fvec=new double[m]; //no need to populate 
        double ftol=1e-08; //tolerance
        double xtol=1e-08; //tolerance
        double gtol=1e-08; //tolerance
        int maxfev=400; //maximum function evaluations
        double epsfcn=1e-08; //tolerance
        double* diag=new double[n]; //some internal thing
        int mode=1; //some internal thing
        double factor=1; // a default recommended value
        int nprint=0; //don't know what it does
        int info=0; //output variable
        int nfev=0; //output variable will store no. of function evals
        double* fjac=new double[m*n]; //output array of jacobian
        int ldfjac=m; //recommended setting
        int* ipvt=new int[n]; //for internal use
        double* qtf=new double[n]; //for internal use
        double* wa1=new double[n]; //for internal use
        double* wa2=new double[n]; //for internal use
        double* wa3=new double[n]; //for internal use
        double* wa4=new double[m]; //for internal use
        QuantLib::MINPACK::lmdif( m, n, x, fvec, ftol, xtol, gtol, maxfev, epsfcn, diag,  mode,  factor, nprint,  &info, &nfev, fjac,ldfjac, ipvt, qtf, wa1, wa2, wa3, wa4);
        //the below is out result. compare it with the values used in fcn function
        double a = x[0]; 
        double b = x[1];
        cout << " mue-Minpack-Fit: " << a << endl;
        cout << " beta-Minpack-Fit: " << b << endl;
        return 0;
    }
    
    int main(int argc, char** argv) 
    {
    
        fitAndTest();
        return 0;
    }
    


  • LOL.

    warum eigentlich der Umweg über Histogramme und nicht direkt maximum likelyhood? Und warum musste es Levenberg-Marquardt sein? Bei dem problem reicht dir auch primitiver Gradeientenabstieg.


Log in to reply