chapter 25 - general linear least squares
tniram
03-10-2010, 01:07 PM
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
davekw7x
03-10-2010, 11:30 PM
...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
MPD78
03-11-2010, 08:43 AM
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
davekw7x
03-11-2010, 07:41 PM
...
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
MPD78
03-12-2010, 07:19 AM
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
tniram
03-12-2010, 09:09 AM
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
davekw7x
03-12-2010, 09:19 AM
...without any errors.
Thank you very much for the update!
Regards,
Dave
davekw7x
03-12-2010, 04:45 PM
...
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
MPD78
03-16-2010, 12:24 PM
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
davekw7x
03-16-2010, 02:40 PM
...
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
MarCn
04-09-2010, 04:50 AM
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.
MPD78
04-27-2010, 09:21 AM
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
davekw7x
04-27-2010, 06:12 PM
...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