HTML automatically generated with rman
Table of Contents
lsq_zero, lsq_accum, lsq_solve - least squares fitting utilities
int lsq_zero (n, mat, vec)
int lsq_accum (n, mat, vec, a, w)
int lsq_solve (n, mat, vec, sol)
int lsq_cfill (n, mat, c, sol)
int n, c;
real mat[n*n], vec[n], sol[n], a[n+1], w;
DescriptionThese routines provide a low level interface to solving linear
least squares problems using the Normal Equations (see e.g. Numerical Recepies,
Press et al.). These routines have the advantage that they don’t need the
extra temporary arrays containing the data, i.e. the Design Matrix, and directly
accumulate the data in the appropriate spots in the symmetric Normal Matrix.
lsq_accum accumulates the data into the normal equations matrix that is
inverted in lsq_solve. Before the data is accumulated, the matrix and
right hand side vector need to be reset using lsq_zero. During accumulation,
a weight, w, can be given to each datapoint (Note: this is a linear weight,
by which the data is multiplied. If you’re used to enter the dispersion,
you should enter 1/sigma^2). On output the matrix mat contains the inverse
of the normal matrix, and hence its diagonal elements the square of the
errors of the fitted parameters. The fitted parameters themselves are returned
in the array sol.
ExampleIn this example a large 2D image matrix is fitted
with an intensity gradient of the form I(x,y)=a+bx+cy:
double mat[3*3],vec[3],a[4],
real data[nx*ny];
lsq_zero(3,mat,vec);
for (ix=0; ix<nx; ix++)
for (iy=0; iy<ny; iy++) {
a[0] = 1.0
a[1] = ix;
a[2] = iy;
a[3] = data[ix*ny+iy];
lsq_accum(3,mat,vec,a,1.0);
}
lsq_solve(3,mat,vec,a);
printf("a+bx+cy: a=%f b=%f c=%f\n", a[0],a[1],a[2]);
lsqfit(3NEMO)
, matinv(3NEMO)
, mpfit(3NEMO)
Peter Teuben
~/src/kernel/misc lsq.c
29-sep-90 created PJT
19-feb-92 updated doc, and properly redfined the weights PJT
Table of Contents