Hi,

I use Fitsvd to fit a polynomial using the NR3 structures of Fitsvd (in fitsvd.h) and fpoly (in fit_examples.h). However, I cannot get it to work, it seems either the Fitsvd constructor or the fitting function has problem. I call the Fitsvd function in the main function as follows:

Fitsvd fitting(x, y, SD, fpoly, TOL);

x,y,SD, and TOL contain values.

The output coefficients are wrong (great numbers) or not assigned.
Do you have ideas of the possible problems and solutions, and any working example?
Thanks.

...Fitsvd... any working example...

//
// Driver for fitsvd
//
// davekw7x

//
// Always include nr3.h first
//
#include "../code/nr3.h"

// Will use a Gaussian random number generator
// to simulate measurement noise. Include
// these to be able to declare a Normaldev object
#include "../code/ran.h"
#include "../code/gamma.h"
#include "../code/deviates.h"

// The following are needed for Fitsvd
#include "../code/svd.h"
#include "../code/fitsvd.h"

//
// A global variable used in the fpoly() function.
//
// Note that you can't simply #include unmodified fit_examples.h
// for this problem since the distribution file sets fpoly_np
// equal to 10, and we will be using a fourth degree polynomial
//
// fpoly_np is set to the number of coefficients of a fourth-degree
// polynomial.
//
Int fpoly_np = 5;
//
// Here's the prototype of fpoly:
//
VecDoub fpoly(const Doub x);


int main(void)
{
Int i;

Int numpoints = 201;
Int npol = fpoly_np;

VecDoub x(numpoints); // Independent variable values
VecDoub y(numpoints); // "Measured" values
VecDoub sig(numpoints); // Estimates of standard deviations of measured values

// A "random" object to generate Gaussian noise is initialized with
// a seed obtained from current time
//
Doub perturb = 0.02; // 2% perturbation of "measured" value
Normaldev normaldev(0.0, perturb, time(0));
//
// For test purposes, might want to use a fixed seed
//Normaldev normaldev(0.0, perturb, 876543321);
//
//

cout << fixed << setprecision(6);
//
// Create a set of data points for a polynomial
//
// The polynomial is 1 + 2x + 3x^2 + 4x^3 + 5x^4
//
// Add some "noise" to the values
//
Doub xstart = 0.0;
Doub xmax = 2.0;
Doub deltax = (xmax - xstart) / (numpoints - 1.0);
Doub xx;
Doub dev;
for (i = 0, xx = xstart; i < numpoints; i++, xx += deltax) {
x[i] = xx;
y[i] = 1.0 + xx * (2.0 + xx * (3.0 + xx * (4.0 + xx * 5.0)));
dev = normaldev.dev();
y[i] += y[i] * dev;
}
//
// For perturb = 0.02, we estimate "measurement error" for all points
// to be 2% of the measured value.
//
for (i = 0; i < numpoints; i++) {
sig[i] = perturb * y[i];
}

//
// Invoke the constructor to bind the data and the function
//
Fitsvd fitsvd(x, y, sig, fpoly);

//
// Calculate the fit
//
fitsvd.fit();


cout << endl << "Polynomial coefficients:" << endl << endl;
cout << setw(11) << "Actual"
<< setw(13) << "Estimated"
<< setw(17) << "Err Estimate"
<< endl;
for (i = 0; i < npol; i++) {
cout << setw(12) << (Doub)(i + 1);
cout << setw(12) << fitsvd.a[i] << " +-";
cout << setw(11) << sqrt(fitsvd.covar[i][i]) << endl;
}
cout << endl << "Number of points = " << numpoints
<< ", Chi-squared = " << fitsvd.chisq << endl;
return 0;
}

//
// Paste the fpoly() function from fit_examples.h here
//


Output from a run on my Linux workstation (GNU g++ version 4.1.2):
Polynomial coefficients:

Actual Estimated Err Estimate
1.000000 0.983243 +- 0.009755
2.000000 2.113134 +- 0.140443
3.000000 2.641075 +- 0.521260
4.000000 4.522804 +- 0.634289
5.000000 4.734538 +- 0.229465

Number of points = 201, Chi-squared = 186.473844


Regards,

Dave