Hello,

I am having some trouble in using the lfit function. I have lets say 8 discrete x,y points with no standard deviations on them and I need to fit a curve to them.

I took the example program and modified it so that I assigned
the values to the vectors x and y

x[0] = 254;
x[1] = 508;
x[2] = 762;
x[3] = 1016;
x[4] = 1269;
x[5] = 1523;
x[6] = 1777;
x[7] = 2031;

y[0] = 46720;
y[1] = 72537;
y[2] = 85964;
y[3] = 93631;
y[4] = 98264;
y[5] = 101299;
y[6] = 103388;
y[7] = 104888;

After this I am not sure as to what must be within the function
void funcs(const DP, Vec_O_dP &)

This function is called by xlfit and by using the function provided in the example for my case I am getting a fit that does not match with what I get from say Excel.

I am using the second edition.

Any help in this regard is appreciated.
Thanks

...some trouble in using the lfit ...
If you are doing linear least squares fitting, you can use the somewhat simpler fit() function rather than lfit().

However...

For polynomials of any order, you can try something like the following example. Just change numTerms to be the number of terms in your polynomial. So for a quadric fit you would use numTerms = 3.

Anyhow...

//
// Example of polynomial least squares approximation
//
// davekw7x
//
#include <iostream>
#include <iomanip>
#include <cmath>
#include "nr.h"

using namespace std;

// Function for least squares polynomial curve fit
void fpoly(const DP x, Vec_O_DP & afunc)
{
int siz = afunc.size();
afunc[0] = 1.0;
for (int i = 1; i < siz; i++) {
afunc[i] = afunc[i-1] * x;
}
}

int main()
{
// Number of points
const int numPoints = 8;

// Number of terms in the approximating polynomial.
// The degree of the polynomialis nTerms-1
const int numTerms = 2;
Vec_DP a(numTerms); // This is where the result is stored

// The following are not used unless you want to test
// goodness of fit.
DP chisq;
Vec_BOOL ia(numTerms);

// Values of sig are kind of irrevelant here, since we
// are not necessarily modeling a polynomial function
// and "goodness of fit" may not mean much.
// The xlfit function will crash if components of sig
// are zero, so set them to something else.
Vec_DP sig(1.0, numPoints);
Mat_DP covar(numTerms,numTerms);


// Define arrays of points that will be used to
// populate the vectors that xlfit needs.
DP xArray[numPoints] = {
254, 508, 762, 1016, 1269, 1523, 1777, 2031
};

DP yArray[numPoints] = {
46720, 72537, 85964, 93631, 98264, 101299, 103388, 104888
};

Vec_DP x(xArray, numPoints);
Vec_DP y(yArray, numPoints);

for (int i = 0; i < numTerms; i++) {
ia[i]=true;
}
NR::lfit(x, y, sig, a, ia, covar, chisq, fpoly);

cout << scientific;
for (int i = 0; i < numTerms; i++) {
cout << " a[" << i << "] = " << a[i];
}
cout << endl << endl;

//
// Use Horners scheme for evaluating the approximating
// polynomial at each value of x.
//
Vec_DP z(numPoints);
for (int i = 0; i < numPoints; i++) {
DP term = a[numTerms - 1];
for (int j = numTerms - 2; j >= 0; j--) {
term = a[j] + term * x[i];
}
z[i] = term;
}

cout << "LS Approx = a[0]";
for (int i = 1; i < numTerms; i++) {
cout << " + a[" << i << "] * x ** " << i;
}
cout << endl << endl;

cout << " X data Y data LS Approx" << endl;
cout << " --------------------------------------------" << endl;
for (int i = 0; i < numPoints; i++) {
cout << setw(15) << x[i]
<< setw(15) << y[i]
<< setw(15) << z[i]
<< endl;
}
return 0;
}


Output using GNU g++ version 4.1.2 on my Centos 5.6 Workstation:

a[0] = 5.553511e+04 a[1] = 2.871008e+01

LS Approx = a[0] + a[1] * x ** 1

X data Y data LS Approx
--------------------------------------------
2.540000e+02 4.672000e+04 6.282747e+04
5.080000e+02 7.253700e+04 7.011983e+04
7.620000e+02 8.596400e+04 7.741219e+04
1.016000e+03 9.363100e+04 8.470455e+04
1.269000e+03 9.826400e+04 9.196820e+04
1.523000e+03 1.012990e+05 9.926056e+04
1.777000e+03 1.033880e+05 1.065529e+05
2.031000e+03 1.048880e+05 1.138453e+05

See attachments for plots of the points and least squares straight line fit and least squares quadratic fit.

Regards,

Dave

Thank you davekw7x for your quick help.