Jochen1980
07-31-2011, 10:48 AM
Hi, I am trying to get Fitmrq to work. The tutorial which is given in the book and other mentioned threads related to Daves example are running. Anyway I am somewhat overstrained how to bend the example to my case:
I want to fit mue and beta of a gumbel distribution to incoming experiment data. I do not expect that all gaps were filled by users here, but I need some more hints to push this snippet in the direction, that it will be running soon.
Thanks in advance.
// user-supplied function of tutorial
void fgauss(const Doub x, VecDoub_I &a, Doub &y, VecDoub_O &dyda) {
Int i, na = a.size();
Doub fac, ex, arg;
y = 0.;
// TODO was passiert da?
for (i = 0; i < na - 1; i += 3)
{
arg = (x - a[i + 1]) / a[i + 2];
ex = exp(-SQR(arg));
fac = a[i] * ex * 2. * arg;
y += a[i] * ex;
dyda[i] = ex;
dyda[i + 1] = fac / a[i + 2];
dyda[i + 2] = fac * arg / a[i + 2];
}
}
// Struktur, welche den Fit-Vorgang verwaltet und steuert
struct Fitmrq {
...
};
void nrCurveFitUserFunction( const Doub x, VecDoub_I &a, Doub &y, VecDoub_O &dyda )
{
// User supplied function for Gumbel, ref en.wikipedia.org/wiki/Gumbel_distribution
// Gumbel-distribution (parameter mue und beta):
// cdf: y = exp( - exp( - ( x - mue ) / beta ) );
// pdf: y = exp( ( -1 * ( x - mue ) ) / beta ) * exp( -1 * ( exp( ( -1 * ( x - mue ) ) / beta )) ) / beta ;
// ---- Tutorial Gauss-Function
// User-provided-function, ref fgauss() in example
Int i, na = a.size();
Doub fac, ex, arg;
y = 0.;
for (i = 0; i < na - 1; i += 3)
{
arg = (x - a[i + 1]) / a[i + 2];
ex = exp(-SQR(arg));
fac = a[i] * ex * 2. * arg;
y += a[i] * ex;
dyda[i] = ex;
dyda[i + 1] = fac / a[i + 2];
dyda[i + 2] = fac * arg / a[i + 2];
}
}
void nrCurveFit() // main
{
// Kapitel 15.5.2; S. 801 ff fitmrq.h
cout << "nrCurveFit()... \n";
// Refs:
// - informative threads:
// numerical-recipes.com/forum/showthread.php?t=1689
// nr.com/forum/archive/index.php/t-2066.html
// nr.com/forum/archive/index.php/t-2085.html
// Model here cumulative Gumbel-distribution
// Gumbel-distribution (parameter mue und beta):
// cdf: y = exp( - exp( - ( x - mue ) / beta ) );
// pdf: y = exp( ( -1 * ( x - mue ) ) / beta ) * exp( -1 * ( exp( ( -1 * ( x - mue ) ) / beta )) ) / beta ;
/*
// data
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);
// Guess
const Doub guessParas[] = {0.2 , 0.002}; // mue, beta
Int numPoints = 30; // number of points
const Int maCf = sizeof (guessParas) / sizeof (guessParas[0]); // I deleted all 'actual'-related expression, because I just got experimental data
VecDoub guessCf(maCf, guessParas);
VecDoub xCf(numPoints);
VecDoub yCf(numPoints);
VecDoub sdCf(numPoints);
Fitmrq mymrqCf(xCf, yCf, sdCf, guessCf, nrCurveFitUserFunction); // bind data
mymrqCf.fit(); // start fit
// print results - form Tutorial
cout << "Initial guess: ";
for (int i = 0; i < guessCf.size(); i++) {
cout << guessParas[i] << " ";
}
cout << endl;
//cout << "noise : " << noise << endl << endl;
cout << scientific << setprecision(5);
cout << setw(10) << "Actual"
<< setw(16) << "Fitted"
<< " "
<< setw(13) << "Err. Est."
<< endl;
cout << fixed;
cout << "Number of points = " << numPoints
<< ", Chi-squared = " << mymrqCf.chisq << endl;
**/
// ---------------
// Tutorial Start
int i;
Doub noise = 0.01; // Standard deviation
Normaldev ran(0, noise, 12345678); // Gaussian noise generator
//Normaldev ran(0, noise, time(0)); // Gaussian noise generator
Int num_points = 101;
const Doub actual_d[] = {5.0, 2.0, 3.0}; // actual guess
const Doub guess_d[] = {4.8, 2.2, 2.8}; // good guess
const Int ma = sizeof (actual_d) / sizeof (actual_d[0]);
VecDoub actual(ma, actual_d);
VecDoub guess(ma, guess_d);
VecDoub x(num_points);
VecDoub y(num_points);
VecDoub sd(num_points);
Doub xstart = 0.0;
Doub xmax = 10.0;
Doub deltax = xmax / (num_points - 1);
Doub xx;
Doub r;
cout << fixed << showpos;
for (i = 0, xx = xstart; i < num_points; i++, xx += deltax) {
x[i] = xx;
r = ran.dev();
y[i] = actual[0] * exp(-SQR((xx - actual[1]) / actual[2]));
y[i] += r * y[i];
sd[i] = noise * y[i];
cout << "y[" << i << "] = " << y[i] << ", r = " << r << ", sd[" << i << "] = " << sd[i] << endl;
}
// -------------------------------------------------------------------------
// Fit vorbereiten
// x = Sequenz der Werte, etwas komisch notiert die for-Schleife;
// y = Verrauschte y-Werte, also Experimentergebnisse;
// sd = in jedem Fall hat der Vektor Punktelaenge, mal schauen
// guess = Startwerte
// fgauss = die user-supplied-function
Fitmrq mymrq(x, y, sd, guess, fgauss); // Daten binden
mymrq.fit(); // Fit starten
// Ergebnis ausgeben und
cout << "Initial guess: ";
for (int i = 0; i < guess.size(); i++) {
cout << guess[i] << " ";
}
cout << endl;
cout << "noise : " << noise << endl << endl;
cout << scientific << setprecision(5);
cout << setw(10) << "Actual"
<< setw(16) << "Fitted"
<< " "
<< setw(13) << "Err. Est."
<< endl;
for (i = 0; i < actual.size(); i++) {
cout << setw(13) << actual[i]
<< setw(16) << mymrq.a[i]
<< " +- " << sqrt(mymrq.covar[i][i])
<< endl;
}
cout << endl;
cout << fixed;
cout << "Number of points = " << num_points
<< ", Chi-squared = " << mymrq.chisq << endl;
//return 0;
// ENDE BEISPIEL
cout << "... // nrCurveFit()\n";
}
I want to fit mue and beta of a gumbel distribution to incoming experiment data. I do not expect that all gaps were filled by users here, but I need some more hints to push this snippet in the direction, that it will be running soon.
Thanks in advance.
// user-supplied function of tutorial
void fgauss(const Doub x, VecDoub_I &a, Doub &y, VecDoub_O &dyda) {
Int i, na = a.size();
Doub fac, ex, arg;
y = 0.;
// TODO was passiert da?
for (i = 0; i < na - 1; i += 3)
{
arg = (x - a[i + 1]) / a[i + 2];
ex = exp(-SQR(arg));
fac = a[i] * ex * 2. * arg;
y += a[i] * ex;
dyda[i] = ex;
dyda[i + 1] = fac / a[i + 2];
dyda[i + 2] = fac * arg / a[i + 2];
}
}
// Struktur, welche den Fit-Vorgang verwaltet und steuert
struct Fitmrq {
...
};
void nrCurveFitUserFunction( const Doub x, VecDoub_I &a, Doub &y, VecDoub_O &dyda )
{
// User supplied function for Gumbel, ref en.wikipedia.org/wiki/Gumbel_distribution
// Gumbel-distribution (parameter mue und beta):
// cdf: y = exp( - exp( - ( x - mue ) / beta ) );
// pdf: y = exp( ( -1 * ( x - mue ) ) / beta ) * exp( -1 * ( exp( ( -1 * ( x - mue ) ) / beta )) ) / beta ;
// ---- Tutorial Gauss-Function
// User-provided-function, ref fgauss() in example
Int i, na = a.size();
Doub fac, ex, arg;
y = 0.;
for (i = 0; i < na - 1; i += 3)
{
arg = (x - a[i + 1]) / a[i + 2];
ex = exp(-SQR(arg));
fac = a[i] * ex * 2. * arg;
y += a[i] * ex;
dyda[i] = ex;
dyda[i + 1] = fac / a[i + 2];
dyda[i + 2] = fac * arg / a[i + 2];
}
}
void nrCurveFit() // main
{
// Kapitel 15.5.2; S. 801 ff fitmrq.h
cout << "nrCurveFit()... \n";
// Refs:
// - informative threads:
// numerical-recipes.com/forum/showthread.php?t=1689
// nr.com/forum/archive/index.php/t-2066.html
// nr.com/forum/archive/index.php/t-2085.html
// Model here cumulative Gumbel-distribution
// Gumbel-distribution (parameter mue und beta):
// cdf: y = exp( - exp( - ( x - mue ) / beta ) );
// pdf: y = exp( ( -1 * ( x - mue ) ) / beta ) * exp( -1 * ( exp( ( -1 * ( x - mue ) ) / beta )) ) / beta ;
/*
// data
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);
// Guess
const Doub guessParas[] = {0.2 , 0.002}; // mue, beta
Int numPoints = 30; // number of points
const Int maCf = sizeof (guessParas) / sizeof (guessParas[0]); // I deleted all 'actual'-related expression, because I just got experimental data
VecDoub guessCf(maCf, guessParas);
VecDoub xCf(numPoints);
VecDoub yCf(numPoints);
VecDoub sdCf(numPoints);
Fitmrq mymrqCf(xCf, yCf, sdCf, guessCf, nrCurveFitUserFunction); // bind data
mymrqCf.fit(); // start fit
// print results - form Tutorial
cout << "Initial guess: ";
for (int i = 0; i < guessCf.size(); i++) {
cout << guessParas[i] << " ";
}
cout << endl;
//cout << "noise : " << noise << endl << endl;
cout << scientific << setprecision(5);
cout << setw(10) << "Actual"
<< setw(16) << "Fitted"
<< " "
<< setw(13) << "Err. Est."
<< endl;
cout << fixed;
cout << "Number of points = " << numPoints
<< ", Chi-squared = " << mymrqCf.chisq << endl;
**/
// ---------------
// Tutorial Start
int i;
Doub noise = 0.01; // Standard deviation
Normaldev ran(0, noise, 12345678); // Gaussian noise generator
//Normaldev ran(0, noise, time(0)); // Gaussian noise generator
Int num_points = 101;
const Doub actual_d[] = {5.0, 2.0, 3.0}; // actual guess
const Doub guess_d[] = {4.8, 2.2, 2.8}; // good guess
const Int ma = sizeof (actual_d) / sizeof (actual_d[0]);
VecDoub actual(ma, actual_d);
VecDoub guess(ma, guess_d);
VecDoub x(num_points);
VecDoub y(num_points);
VecDoub sd(num_points);
Doub xstart = 0.0;
Doub xmax = 10.0;
Doub deltax = xmax / (num_points - 1);
Doub xx;
Doub r;
cout << fixed << showpos;
for (i = 0, xx = xstart; i < num_points; i++, xx += deltax) {
x[i] = xx;
r = ran.dev();
y[i] = actual[0] * exp(-SQR((xx - actual[1]) / actual[2]));
y[i] += r * y[i];
sd[i] = noise * y[i];
cout << "y[" << i << "] = " << y[i] << ", r = " << r << ", sd[" << i << "] = " << sd[i] << endl;
}
// -------------------------------------------------------------------------
// Fit vorbereiten
// x = Sequenz der Werte, etwas komisch notiert die for-Schleife;
// y = Verrauschte y-Werte, also Experimentergebnisse;
// sd = in jedem Fall hat der Vektor Punktelaenge, mal schauen
// guess = Startwerte
// fgauss = die user-supplied-function
Fitmrq mymrq(x, y, sd, guess, fgauss); // Daten binden
mymrq.fit(); // Fit starten
// Ergebnis ausgeben und
cout << "Initial guess: ";
for (int i = 0; i < guess.size(); i++) {
cout << guess[i] << " ";
}
cout << endl;
cout << "noise : " << noise << endl << endl;
cout << scientific << setprecision(5);
cout << setw(10) << "Actual"
<< setw(16) << "Fitted"
<< " "
<< setw(13) << "Err. Est."
<< endl;
for (i = 0; i < actual.size(); i++) {
cout << setw(13) << actual[i]
<< setw(16) << mymrq.a[i]
<< " +- " << sqrt(mymrq.covar[i][i])
<< endl;
}
cout << endl;
cout << fixed;
cout << "Number of points = " << num_points
<< ", Chi-squared = " << mymrq.chisq << endl;
//return 0;
// ENDE BEISPIEL
cout << "... // nrCurveFit()\n";
}