Hello!

Who knows what sig refers to in fitsvd.h or fitlin.h?
Could you give an example?
Is there any examples for showing how to use nr3?
Thanks!
xi shi.

Who knows what sig refers to in fitsvd.h or fitlin.h?
As mentioned in the early sections of chapter 15, if you have an array of "measurement error" values for the data points and if the "measurement errors" have a gaussian distribution, the covariance matrix can give an idea of the accuracy of the calculated parameters, and the value of chi-square can give a sense of "goodness of fit." The elements of the array sig are used to feed your estimates of measurement errors to the fitsvd constructor. (The identifier sig refers to the standard deviation of the error value, or sigma.)

If you don't have a reasonable guess of the magnitudes of measurement errors, you can just use values proportional to the values of the data points and you still may be able to get a reasonable set of parameters (but you will probably lose that warm fuzzy feeling that you might have seen about the goodness of fit).


Could you give an example?Using a fourth degree polynomial, as mentioned in section 15.4.3 with the svd fitting routine of 15.4.2:


//
// Driver for test of fitsvd
// davekw7x
// May, 2008
//
#include "../code/nr3.h"
#include "../code/ran.h"
#include "../code/gamma.h"
#include "../code/deviates.h"
#include "../code/svd.h"
#include "../code/fitsvd.h"

/* Function to evaluate 4th degree polynomial */
VecDoub fpoly4(const Doub x)
{
Int j;
VecDoub p(5);
p[0] = 1.0;
for (j = 1; j < 5; j++)
p[j] = p[j - 1] * x;
return p;
}


int main(void)
{
const int NPT = 101; // Number of points of sample data
const int NPOL = 5; // Fourth degree polynomial has five coefficients
//
// Will construct an array of points perturbed by gaussian noise
// that simulates measurement errors.
//
const Doub NOISE = 0.002; // Noise standard deviation will be 0.2% of measured value
int i;

VecDoub x(NPT), y(NPT), sig(NPT), a(NPOL), w(NPOL);
MatDoub cvm(NPOL, NPOL), u(NPT, NPOL), v(NPOL, NPOL);

//
// Declare a random object.
//
// To get a different noise each run, seed with time(0).
//
// For test and comparison purposes it is sometimes useful
// to use the same seed for all runs.
//

//Normaldev gauss(0.0, NOISE, time(0));
Normaldev gauss(0.0, NOISE, 123456);

//
// Generate samples, perturbed by gausian noise
//
// The polynomial is 1 + 2x + 3x^2 + 4x^3 + 5x^4
//
//
//We will have NPT values of x on [0,2]
//
Doub deltax = 2.0/(NPT-1);
Doub xvalue = 0;

for (i = 0; i < NPT; xvalue+=deltax, i++) {
x[i] = xvalue;
y[i] = 1.0 + x[i] * (2.0 + x[i] * (3.0 + x[i] * (4.0 + x[i] * 5.0)));
y[i] += gauss.dev()*y[i];
}
//
// Use same standard error for all sig[i] values.
// Start with a reasonable estimate, based on knowledge
// of "measurement error" noise
//
// Try it again with NOISE_FACTOR = 1 and see what happens.
//
double NOISE_FACTOR = NOISE;
for (i = 0; i < NPT; i++) {
sig[i] = NOISE_FACTOR*y[i];
}

//
// Create the object
//
Fitsvd myFitsvd(x, y, sig, fpoly4);

//
// Do the fit
//
myFitsvd.fit();

// If the estimate of measurement errors was "good", and
// the data are actually modeled reasonably by the
// function, the covariance matrix can give a good idea of
// the accuracy of each coefficient.
//
// Furthermore, we expect to see that the value of chi-square
// is on the order of NPT-NPOL (which is the number of degrees
// of freedom).

cout << fixed << setprecision(6);
cout << endl << "fitsvd for 1 + 2x + 3x^2 + 4x^3 + 5x^4:" << endl << endl;
cout << "With noise factor = " << NOISE_FACTOR << endl;

cout << setw(12) << "Coefficient"
<< setw(10) << "Approx"
<< setw(22) << "sqrt(covariance)" << endl;

for (i = 0; i < NPOL; i++) {
cout << setw(10) << static_cast < double >(i + 1);
cout << setw(13) << myFitsvd.a[i] << " +-";
cout << setw(12) << sqrt(myFitsvd.covar[i][i]) << endl;
}
cout << endl << "Chi-squared " << setw(12) << myFitsvd.chisq << endl;

return 0;
}


Output using GNU g++ version 4.1.2 on my Linux system:


fitsvd for 1 + 2x + 3x^2 + 4x^3 + 5x^4:

With noise factor = 0.002000
Coefficient Approx sqrt(covariance)
1.000000 1.000739 +- 0.001307
2.000000 1.990259 +- 0.019124
3.000000 3.021520 +- 0.071920
4.000000 4.001873 +- 0.088166
5.000000 4.987233 +- 0.032030

Chi-squared 107.155592


This used an array of error values that were calculated from the standard deviation of the measurement noise.

A cursory examination shows:

The coefficients look reasonable, the error estimates are pretty good, and the value of chi-square is on the order of the number of degrees of freedom of the system (101 points, 5 polynomial coefficients give 96 degrees of freedom).

Now, If I change it so that sig[i] = y[i], here's what I see:

fitsvd for 1 + 2x + 3x^2 + 4x^3 + 5x^4:

With noise factor = 1.000000
Coefficient Approx sqrt(covariance)
1.000000 1.000739 +- 0.653286
2.000000 1.990259 +- 9.562206
3.000000 3.021520 +- 35.959862
4.000000 4.001873 +- 44.083007
5.000000 4.987233 +- 16.015072

Chi-squared 0.000429


Since I know what polynomial was to use to generate the data points, I can see at a glance that the coefficients are OK.

However...
The other things no longer give me any idea at all of goodness of fit. If the measurement errors weren't gaussian, I may or may not have obtained reasonable coefficient approximations, but I would definitely need some other way of guessing at goodness of fit.

Regards,

Dave

Thank you, Dave