General Linear Least Squares


LoneJoe
03-14-2007, 04:18 PM
Hi,

I am trying to use the methods described in 15.4 General Least Squares, so create a polynomial best fit equation based on X and Y data.

Keeping it simple, I am using only four data points, and checking the answers with Excel. For example:
X, Y
1, -0.7
2, -8.6
3, -20.9
4, -35.8

The answer should be y = 0.3x^3 - 4x^2 + 2x +1; and therefore I would expect my array pointed to by the 'a' function in the svdfit function should return 0.3, -4, 2, and 1.

Well, beyond my expections, I am not getting anywhere.

The function requires a 'basis function'. What is that?

Is there a different function I should be calling? For example, I can understand defining the matrix scratch pad areas (as opposed to having them dynamically created in the function); but is there something else I should do?

Does anyone have an simple X,Y example that I can see and build from?

LoneJoe
03-15-2007, 10:12 AM
I found the example code for svdfit.

How can I modify it for a simple example?
For example, with the following set of data:
X, Y
1 4.5
2 9
3 17.5
4 30
5 46.5


I tried the following,.. but the a values returned are no where near what I expected (y = 2x2 - 1.5x + 4) as I get a bunch of zeros... what am I doing wrong?

/* Driver for routine svdfit */

#include <stdio.h>
#include <math.h>
#define NRANSI
#include "nr.h"
#include "nrutil.h"

#define NPT 4
#define SPREAD 0.02
#define NPOL 5

int main(void)
{
long idum=(-911);
int i;
float chisq,*x,*y,*sig,*a,*w,**cvm,**u,**v;

x=vector(1,NPT);
y=vector(1,NPT);
sig=vector(1,NPT);
a=vector(1,NPOL);
w=vector(1,NPOL);
cvm=matrix(1,NPOL,1,NPOL);
u=matrix(1,NPT,1,NPOL);
v=matrix(1,NPOL,1,NPOL);
/*
for (i=1;i<=NPT;i++) {
x[i]=0.02*i;
y[i]=1.0+x[i]*(2.0+x[i]*(3.0+x[i]*(4.0+x[i]*5.0)));
y[i] *= (1.0+SPREAD*gasdev(&idum));
sig[i]=y[i]*SPREAD;
} */
x[0] = 1;
x[1] = 2;
x[2] = 3;
x[3] = 4;
x[5] = 5;

sig[0] = .100;
sig[1] = .100;
sig[2] = .100;
sig[3] = .100;
sig[4] = .100;

x[0] = 4.5;
x[1] = 9;
x[2] = 17.5;
x[3] = 30;
x[4] = 46.5;

svdfit(x,y,sig,NPT,a,NPOL,u,v,w,&chisq,fpoly);
svdvar(v,NPOL,w,cvm);
printf("\npolynomial fit:\n\n");
for (i=1;i<=NPOL;i++)
printf("%12.6f %s %10.6f\n",a[i]," +-",sqrt(cvm[i][i]));
printf("\nChi-squared %12.6f\n",chisq);
svdfit(x,y,sig,NPT,a,NPOL,u,v,w,&chisq,fleg);
svdvar(v,NPOL,w,cvm);
printf("\nLegendre polynomial fit:\n\n");
for (i=1;i<=NPOL;i++)
printf("%12.6f %s %10.6f\n",a[i]," +-",sqrt(cvm[i][i]));
printf("\nChi-squared %12.6f\n",chisq);
free_matrix(v,1,NPOL,1,NPOL);
free_matrix(u,1,NPT,1,NPOL);
free_matrix(cvm,1,NPOL,1,NPOL);
free_vector(w,1,NPOL);
free_vector(a,1,NPOL);
free_vector(sig,1,NPT);
free_vector(y,1,NPT);
free_vector(x,1,NPT);
return 0;
}