Regarding using a vector of standard deviations in a linear fit using SVD, chapter 15 says that by defining the model matrix values as

A_ij = x_ij / sig_i ( as in 15.4.4, but I'm not using basis functions here )

and the response values as

b_i = y_i / sig_i ( as in 15.4.5 )

where "y" is the response value vector, and "sig" the vector containing their estimated standard errors.

Using these as input data to svd I should be able to get the best-fit parameter values, according to eq. 15.14.16


For the following simple test I get the following:

I define a simple two-parameter linear model an set their values to the following ones: [1, 0.3]
I define a model matrix as
[1,0]
[0,1]
[0,1]
so that its corresponding response is
y = [ 1, 0.3, 0.3 ]

after svd of the model matrix, I solve for "y" and get exactly the starting parameters [1, 0.3], as expected ( I'm using a translation of svd to a scripting language, not Fitsvd itself)


Then I decide to introduce dummy standard deviations for such response values. If I define

sig = [1, 1, 1]

and modify the response and model matrix as indicated above to get "A" and "b", and I get exactly the same results as before, as expected.

Then, I try other "sig" vectors

For sig = [0.1, 0.1, 0.1] , I get [10, 3] as best fit parameters.

For sig = [1, 1, 0.1] , I get [1, 1.65] as best fit parameters.

For sig = [0.1, 0.1, 0.01] , I get [10, 16.5] as best fit parameters.

I suspect I'm incorrectly interpreting sig values as being used as weights and such thing is confusing me (assuming normally distributed errors, it would be unlikely that measurements widely differing in standard deviation would yield the same response)

Do these results make sense or I have a bug in my code?
Any better example to try?
Some hint on what's going on here?
I'd like to use this method on data where the "sig" values of the "y" values can be measured empirically.

- many thanks in advance

...Do these results make sense...
Not to me:


//
// Demo of use of fitsvd for linear regression
//
// davekw7x March, 2010

#include "../code/nr3.h"
#include "../code/svd.h"
#include "../code/fitsvd.h"


//
// Basis for
// y = b0 * x[0] + b1 * x[1] + ...
//
VecDoub funks(VecDoub_I & x)
{
VecDoub p(x);
return p;
}


int main()
{
//
// The left-hand sides
//
Doub x_array[3*2] = {
1.0, 0.0,
0.0, 1.0,
0.0, 1.0
};
//
// The right-hand sides
//
Doub y_array[3] = {
1.0,
0.3,
0.3
};
//
// Try several sets of ssig values
//
Doub sig_array[5][3] = {
{1.0e-50, 1.0e-50, 1.0e-50}, // Can't use all zeros, so use small
{ 1.0, 1.0, 1.0},
{ 0.1, 0.1, 0.1},
{ 1.0, 1.0, 0.1},
{ 0.1, 0.1, 0.01}
};

MatDoub x(3,2, x_array);
VecDoub y(3, y_array);

cout << "------------------------------------------------"
<< endl << endl;


for (int k = 0; k < 5; k++) {
VecDoub ssig(3, sig_array[k]);

Fitsvd fitsvd(x, y, ssig, funks); // Instantiate the object
fitsvd.fit(); // Do the work

cout << scientific << setprecision(5);
cout << "ssig: ";
for (int i = 0; i < ssig.size(); i++) {
cout << " " << ssig[i];
}
cout << endl << endl;
for (int i = 0; i < fitsvd.a.size(); i++) {
cout << " b" << i << " = ";
cout << setw(13) << fitsvd.a[i] << " +-";
cout << setw(13) << sqrt(fitsvd.covar[i][i]) << endl;
}
cout << endl;
cout << "------------------------------------------------"
<< endl << endl;
}

return 0;
}


Output (GNU g++ version 4.1.2):

------------------------------------------------

ssig: 1.00000e-50 1.00000e-50 1.00000e-50

b0 = 1.00000e+00 +- 1.00000e-50
b1 = 3.00000e-01 +- 7.07107e-51

------------------------------------------------

ssig: 1.00000e+00 1.00000e+00 1.00000e+00

b0 = 1.00000e+00 +- 1.00000e+00
b1 = 3.00000e-01 +- 7.07107e-01

------------------------------------------------

ssig: 1.00000e-01 1.00000e-01 1.00000e-01

b0 = 1.00000e+00 +- 1.00000e-01
b1 = 3.00000e-01 +- 7.07107e-02

------------------------------------------------

ssig: 1.00000e+00 1.00000e+00 1.00000e-01

b0 = 1.00000e+00 +- 1.00000e+00
b1 = 3.00000e-01 +- 9.95037e-02

------------------------------------------------

ssig: 1.00000e-01 1.00000e-01 1.00000e-02

b0 = 1.00000e+00 +- 1.00000e-01
b1 = 3.00000e-01 +- 9.95037e-03

------------------------------------------------


Regards,

Dave

davekw7x,

I am trying to duplicate the results you posted but I keep getting an access violation error from svd.h.

I verified that my svd.h routine is running correctly for a square matrix following your example posted here.

http://www.nr.com/forum/showthread.php?t=1437&highlight=svd

The access violation error occurs in this function from svd.

void SVD::solve(VecDoub_I &b, VecDoub_O &x, Doub thresh = -1.0)

at the loop that performs the matrix multiplication by V to get answer.

for(jj=0;jj<n;j++) s += v[j][jj]*tmp[jj];

I am not sure what I am doing wrong. I copied and pasted your example code to avoid typos but I still get the access violation error.

I can't find any errors that I might have caused in fitsvd.h.

If you can point me towards destiny, I would be very happy.

Thanks
Matt

...
I am not sure what I am doing wrong. I copied and pasted your example code to avoid typos but I still get the access violation error....I am using "virgin" fitsvd.h from the NR3 CD. It is exactly the same as is in the first printing of the NR3 text (As Far As I Know).

The test output that I showsed is the result of using GNU g++ version 4.1.2 on a 32-bit Centos 5.4 Linux workstation.

Is your code hand-entered (copied from the text) or NR software download or NR3 CD or what? What compiler and operating system are you using?


Regards,

Dave

I am using the code from the electronic edition of the NR3 book.

It is hand typed in by me from the electronic edition. I have checked for typos but there aren't any. (I know ..... )

I am using the 2008 C++ Express edition

Microsoft Visual Studio 2008
Version 9.0.30729.1 SP
Microsoft .NET Framework
Version 3.5 SP1

The error must lie in the fitsvd.h routine. I arrive at this conclusion because I can use svd.h and duplicate your results from the above linked driver example that you provided for another user.

Thanks
Matt

Not to me:

Dave

Thanks for your time, your piece of advice was very useful to me.

I managed to find the bug in my code causing all the nonsense. Now I can reproduce your results exactly using my code and using virgin NR3 code too, either through Fitsvd.h or using svd.h directly, without any errors.


Best regards

...without any errors.

Thank you very much for the update!

Regards,

Dave

...
I am using the 2008 C++ Express editionWhen I finally got around to testing it with the 2008 version of Visual C++ Express on my Windows XP workstation, I found that the results were the same (except, of course, that the exponent values in the output showed three decimal digits instead of two.)

The point is that I always like to check these things on different compilers because it has happened that sometimes I get different results.

In this case, I guess it's the code that is the problem, as you have already realized.

Regards,

Dave

davekw7x,

I tried the other example that you posted here.

http://www.nr.com/forum/showthread.php?t=1966 (post #6)

I got the same error as the other example program.

Is there something that I have do with svd.h to handle a rectangular matrix? The reason I ask that question is because svd.h works fine for any (all that I have tried) square matrices.

Thanks
Matt

...
Is there something that I have do with svd.h to handle a rectangular matrix?...
In my previous post I mentioned that I am using "virgin" fitsvd.h

However, svd.h had been diddled a little, and I didn't realize that it might be the problem.

In particular, I note that Fitsvd::fit() calls SVD::solve().

Have you made the changes to svd.h that were prescribed in this thread? http://www.nr.com/forum/showthread.php?t=1358

If not, maybe that's the problem.

Regards,

Dave

Thanks for your time, your piece of advice was very useful to me.

I managed to find the bug in my code causing all the nonsense. Now I can reproduce your results exactly using my code and using virgin NR3 code too, either through Fitsvd.h or using svd.h directly, without any errors.


Best regards

Thanks very much for the update, really helped me get some things straighen out.

Dave,

I have been side tracked lately and I am just now getting back to this issue.

Have you made the changes to svd.h that were prescribed in this thread? http://www.nr.com/forum/showthread.php?t=1358

Yes, I have this code in the svd.h routine but I still cannot obtain results from svd.h for a rectangular matrix. Any other thoughts as to what might be causing this problem?

Thanks
Matt

...what might be causing this problem?My guess: Bad copy of code.

Here's a test program that uses svd.h and shows all output matrices:

#include "../code/nr3.h"
#include "../code/svd.h"

//
// Driver for routine svdcmp
//

int main()
{
Int i, j, k;

const Int m = 4; // Number of rows
const Int n = 3; // Number of columns;
Doub a_array[m * n] = { // Array used to initialize the MatDoub
4.0, 3.0, 5.0,
2.0, 5.0, 8.0,
3.0, 6.0, 10.0,
4.0, 5.0, 11.0
};

MatDoub a(m, n, a_array);
VecDoub w(n);
MatDoub u(m, n);
MatDoub v(n, n);


SVD svd(a);
cout << "*******After instantiation of svd*******" << endl << endl;
cout << "Matrix U (svd.u)" << endl;

cout << fixed << setprecision(6);
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << svd.u[i][j];
}
cout << endl;
}
cout << endl;
cout << "Diagonal of matrix W (svd.w)" << endl;
for (i = 0; i < n; i++) {
cout << setw(12) << svd.w[i];
}
cout << endl << endl;
cout << "Matrix V-transpose (svd.v)" << endl;
for (i = 0; i < n; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << svd.v[i][j];
}
cout << endl;
}
cout << endl;
cout << "Check product against original matrix:" << endl << endl;

cout << "Product U*W*(V-transpose):" << endl;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
a[i][j] = 0.0;
for (k = 0; k < n; k++) {
a[i][j] += svd.u[i][k] * svd.w[k] * svd.v[j][k];
}
cout << setw(12) << a[i][j];
}
cout << endl;
}
cout << endl;
cout << "Original matrix:" << endl;
for (i = 0; i < m; i++) {
for (j = 0; j < n; j++) {
cout << setw(12) << a[i][j];
}
cout << endl;
}
cout << endl;
return 0;
}


Output

*******After instantiation of svd*******

Matrix U (svd.u)
0.321103 -0.855318 0.397029
0.455375 0.411113 0.362370
0.570888 0.288215 0.314183
0.603004 -0.127866 -0.782524

Diagonal of matrix W (svd.w)
21.049381 2.370210 1.142656

Matrix V-transpose (svd.v)
0.300239 -0.947539 0.109666
0.459896 0.244522 0.853642
0.835675 0.205861 -0.509185

Check product against original matrix:

Product U*W*(V-transpose):
4.000000 3.000000 5.000000
2.000000 5.000000 8.000000
3.000000 6.000000 10.000000
4.000000 5.000000 11.000000

Original matrix:
4.000000 3.000000 5.000000
2.000000 5.000000 8.000000
3.000000 6.000000 10.000000
4.000000 5.000000 11.000000



Regards,

Dave