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
-
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, ¤tRealClassLower, ¤tRealClassUpper ); 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, ¤tRealClassLowerCumulated, ¤tRealClassUpperCumulated ); 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.
