Re: Polynomial Fit Program

Ray Koopman wrote:
On Dec 6, 9:56 am, Helmut Jarausch <jarau...@xxxxxxxxx> wrote:
Thomas Smid wrote:
Hi,
Does anybody know an open source code for a one dimensional polynomial
(n-th order) least squares fit program (preferably C++ or Fortran)?
Something similar to the Java applet under
http://www3.sympatico.ca/mcomeau/webpublic/javapage/reg/reg.htmwould
do, but I need the source code, as I need to modify it to perform a
weighted fit (if this isn't built in already). Alternatively, if
somebody knows the contact details of the author of this particular
Java applet (Michel Comeau), it might help me as well.
I don't know if it helps, but free "Matrix Packages" like Octave or Scilab
have it built in.

What you really need is a Q-R-decomposition from any linear algebra package in C++.

In the simplest case (polynomials of moderate degree), if your data is contained
in the two column vectors x and y, build the matrix

A=[ones(x),x,x.^2,x.^3] // for a polynomial of degree 3
This is Matlab/Scilab/Octave jargon :

ones(x) : a column vector of as many ones as x has components

x.^n : (note the dot in front of ^) a column vector which has
x_i^n as i-th component

Now solve the least squares problem A*c ~ y and c contains the coefficients
of the polynomial (lowest first).

If you have weights, e.g. w_i, then first formally
build the diagonal matrix W= diag(w_1,...) and
solve the least squares problem W*A*c ~ W*y

In practice this is done by multiplying the i-th row of A and the
i-th component of y by w_i.

As you can see, if you have a standard (i.e. non weighted) least squares
solver you don't need to modify it.

--
Helmut Jarausch

Lehrstuhl fuer Numerische Mathematik
RWTH - Aachen University
D 52056 Aachen, Germany

Solving W*A*c ~ W*y will weight (A*c - y).^2 by W^2.
Is that what you take "weighted by W" to mean?

Yes,

given the i-th component of the residual r_i:= e_i^T *(A*c-y)
this minimizes sum_i { (w_i*r_i)^2 }

Assuming the residuals have expectation 0 and variance sigma_i^2, taking
w_i= 1/sigma_i is the right choice.



--
Helmut Jarausch

Lehrstuhl fuer Numerische Mathematik
RWTH - Aachen University
D 52056 Aachen, Germany
.

Relevant Pages

  • On limitations of the Mautsch Principle (was Re: ineqality)
    ... > squares. ... I would prove the homogeneous polynomial inequality ... polynomials (which would necessarily have to be linear combinations ... So symmetry on the way in does not imply symmetry on the way out. ...
    (sci.math)
  • Re: Polynomial Fit Program
    ... Does anybody know an open source code for a one dimensional polynomial ... least squares fit program? ... In the simplest case (polynomials of moderate degree), ...
    (sci.math.num-analysis)
  • Proposed CPAN module Statistics::LineFit
    ... statistics, the predicted y values and the residuals of the y ... weights are input in a separate array. ... method to check the status of the regression. ... regress- Do the least squares line fit ...
    (comp.lang.perl.modules)
  • Re: Expressing a polynomial as a sum of squares
    ... > sum of squares of polynomials when such a representation exists? ...
    (sci.math.research)
  • Re: Polynomial Fit Program
    ... Does anybody know an open source code for a one dimensional polynomial ... least squares fit program? ... In the simplest case (polynomials of moderate degree), ...
    (sci.math.num-analysis)
Loading