Dear All

I am in search of a "working version of svdfit() in C++". I am having ample of difficulty in executing this function. In case anybody is ready to share the code or experiences would be highly appreciated.

thanking you in advanced!!

Regards
Roy

...working version of svdfit() in C++...

Is your question about some aspect of the function and/or its implementation? Or is it about creating a test program that uses the function with a particular data set? Or about interpreting/using some results from the function? Or what?

In other words: What, exactly, would you like to know?

Regards,

Dave

Dear Dave

At first I thank you for your reply, secondly my basic aim is to have a multidimentional (fits)regression. The way I am trying to solve it, is by replacing the "xi" by vector "xi" in the General Linear Least Squares method which is mentioned in the book.

I tried to do it, but somehow it doesnt work. So my requirement for now is to have a test program that uses the function with a particular data set (A multidimensional regression test program would be the best).

To your notice I am a average mathematician and programmer.

Regards
Roy

...in the General Linear Least Squares method which is mentioned in the book...


#include <iostream>
#include <iomanip>
#include "nr.h"

int NX = 10; // Number of x values
int NY = 10; // Number of y values
int NPT = NX*NY; // Total number of points

//
// For printing out coefficient terms
//
const int NTERMS = 6; // Number of terms in quadratic
char *label[6] = {
" ",
" * x ",
" * y ",
" * x*x",
" * x*y",
" * y*y"
};

//
// xx will be a matrix with NPT rows: Each row contains x,y for a point
//
Mat_DP xx(NPT,2);


//
// Use the "fake" value to index into the global matrix xx
// To calculate the value at each point
//
void quadratic2d(const DP fake, Vec_O_DP & ans)
{
DP x = xx[int(fake)][0];
DP y = xx[int(fake)][1];
ans[0] = 1;
ans[1] = x;
ans[2] = y;
ans[3] = x*x;
ans[4] = x*y;
ans[5] = y*y;
}

int main()
{
// Generate points over this region
DP xmax = 1.0;
DP xmin = 0.0;
DP ymax = 1.0;
DP ymin = 0.0;


DP deltax = (xmax-xmin)/(NX-1.0);
DP deltay = (ymax-ymin)/(NY-1.0);
DP x;
DP y;

const DP NOISE = 0.02; // 2% "measurement noise"

// For test purposes, may want to use same value each time
int idum=(-54321);

// In general, use a different noise set each time, so use the following
//int idum=(-time(0));

int i, j;

DP chisq;

Vec_DP z(NPT); // The function value at each point
Vec_DP sig(NPT); // Measurement errors

Vec_DP a(NTERMS); // Output of svdfit: the coefficients;
Mat_DP u(NPT,NTERMS); // Result from svd decomposition
Mat_DP v(NTERMS,NTERMS); // Result from svd decomposition
Vec_DP w(NTERMS); // Result from svd decomposition
Mat_DP cvm(NTERMS,NTERMS); // Will hold the covariance matrix

// Allow multidimensional fit, as mentioned at end of section 15.4
Vec_DP fakeit(NPT);
for (i = 0; i < NPT; i++) {
fakeit[i] = i;
}

//
// Generate points for z = x^2 + y^2
//
x = xmin;
int rowno = 0;
for (i = 0; i < NX; i++) {
y = ymin;
for (j = 0; j < NY; j++) {
xx[rowno][0] = x;
xx[rowno][1] = y;
z[rowno] = x*x+y*y;
if (!z[rowno]) { // Note that 1/sig[i] is used in svdfit
sig[rowno] = NOISE;
}
else {
sig[rowno] = z[rowno] * NOISE;
}
z[rowno] *= (1.0 + NOISE * NR::gasdev(idum));
y += deltay;
++rowno;
}
x += deltax;
}

// Do the fit
NR::svdfit(fakeit, z, sig, a, u, v, w, chisq, quadratic2d);

// Calculate covariances
NR::svdvar(v, w, cvm);
cout << fixed << setprecision(6);
cout << "Quadratic fit for " << NPT << " points:" << endl << endl;
cout << " Coeff Err est Term" << endl << endl;
for (i = 0; i < NTERMS; i++) {
cout << "("
<< setw(9) << a[i] << " +-"
<< setw(9) << sqrt(cvm[i][i])
<< ")"
<< label[i] << endl;
}
cout << endl << "Chi-squared " << setw(12) << chisq << endl;
return 0;
}


This uses points from z = x^2 + y^2, perturbed by Gaussian "measurement noise," and fits to a generic quadratic of the form

a0 + a1*x + a2*y + a3*x^2 + a4*x*y + a5*y^2


Output, showing error estimate of each term:


Quadratic fit for 100 points:

Coeff Err est Term

( 0.000682 +- 0.000397)
(-0.008752 +- 0.003481) * x
(-0.005100 +- 0.003481) * y
( 1.007710 +- 0.007061) * x*x
( 0.008876 +- 0.011766) * x*y
( 1.006488 +- 0.007061) * y*y

Chi-squared 82.912970


Regards,

Dave

Footnote: It's a lot less ugly with the Fitsvd class in Numerical Recipes Version 3. (No global matrix or funny argument to index into it.)

Hello Dave

Excellent job!! I was waiting for this since long. I am also working on same issue of multidimensional fit. Unfortunately I have all header files and the associated functions which you have used in this program in C, so I am afraid that I may not be able to use your program. Can you or Roy please provide me the C++ version of "nr.h" and "nrutil.h" header file?

I think I can manage the functions :D

thanks!!

Cheers
Niko Näther

Can you or Roy please provide me the C++ version of "nr.h" and "nrutil.h" header file?


The file nr.h is a copyrighted and is included in the Numerical Recipes distribution, along with everything else. The file nrutil.h has been placed in the public domain, and can be downloaded from here: http://www.nr.com/public-domain.html. You will need nrtypes.h also. (Download nrtypes_nr.h and nrutil_nr.h and rename them if you want.)

There is a discussion here:http://www.numerical-recipes.com/forum/showthread.php?t=42

If you have the actual functions from the book, you can compile my example program with the following changes at the beginning:


#include <iostream>
#include <iomanip>

//#include "nr.h"
// nr.h includes the following public domain files, which can
// be obtained from http://www.nr.com/public-domain.html
// (Download nrtypes_nr.h and nrutil_nr.h and rename them.)
#include "nrtypes.h"
#include "nrutil.h"

// nr.h also contains prototypes of functions in the NR namespace
//
// For this program, use the following NR:: function declarations
//
namespace NR
{
DP gasdev(int &idum);
void svdfit(Vec_I_DP &x, Vec_I_DP &y, Vec_I_DP &sig, Vec_O_DP &a,
Mat_O_DP &u, Mat_O_DP &v, Vec_O_DP &w, DP &chisq,
void funcs(const DP, Vec_O_DP &));
void svdvar(Mat_I_DP &v, Vec_I_DP &w, Mat_O_DP &cvm);
//
// Other function prototypes (and their implementations) will be needed for
// functions called from these.
//
};
.
.
.



Regards,

Dave

Footnote: The Numerical Recipes function svdfit() calls svdcmp() and svbksb() so you will need them also. There are a few others, like pythag() and maybe more. Bottom line: If you don't have the Numerical Recipes distribution, there will be a fair amount of work to get everything to compile and execute. It could be done by typing in stuff as needed from the book, but if at all possible, I think you might find it better to get an official distribution.

Thaks a lot Dave

I will try to excute the program, in case if I will have any question. I will contact you.

Regards

Dear Dave

I bought the C++ source code edition and was pretty succesful in executing the program provided by you.
Recently I had an idea of generalising the program which you have provided. My idea is to get a program which does regression, in a sense that it should have “n” number of independent variables, with “p” number power (fixed functions), and nbr number of data points. And then it should ask at the start of the execution, for the same parameters. So in case I have 4 independent variables and 2 fixed function or I had 5 independent variables and 3 fixed function.. Both can be executed. in a single routine as per the data given by the user.

Also I am planning to give the data points through an array to the main routine.

Can anybody please figure it out??

Thanks

Hi

It is interesting for my research project, to know how can we carry on generalisation for linear multivariate regression. Can anybody share me some ideas?

Merci

Scientificgal