Discussion:
better polyfit function
Todd Elliott
2007-09-05 05:40:29 UTC
Permalink
Hope this is the correct place...

Please find a replacement for the current polyfit in Octave attached. I
believe it produces coefficients that are orders of magnitude better than
the current Octave polyfit function. I was told to submit it to the
mailing lists as a bug, hopefully this is the right place.

Thanks
-Todd
John W. Eaton
2007-09-13 19:13:55 UTC
Permalink
On 4-Sep-2007, Todd Elliott wrote:

| Hope this is the correct place...
|
| Please find a replacement for the current polyfit in Octave attached. I
| believe it produces coefficients that are orders of magnitude better than
| the current Octave polyfit function. I was told to submit it to the
| mailing lists as a bug, hopefully this is the right place.
|
| Thanks
| -Todd

What is the origin of this code? Is it derived from some other
source? Also, is there some way to avoid the for loops and vectorize
it?

Before it could be included in Octave as a replacement for the current
polyfit function, it would need to better match the coding conventions
of the rest of Octave, have a doc string (could use the one from the
existing polyfit function), copyright statement, etc.

Also, I think we would want to some examples that show that it
produces correct results (tests would be nice) and some performance
numbers.

jwe
Todd Elliott
2007-09-13 22:06:12 UTC
Permalink
Hi John,
I completely agree with your concerns. As I think I mentioned in my
original e-mail to you, I aquired this code in the Java language some
years ago. I believe it was "freeware" as in no-license-whatsoever from
the following (now defunct link)
<http://www.intergalact.com/java/polyfit/polyfit.html><http://www.intergalact.com/java/polyfit/polyfit.html>
I have used the code in Java and C forms for many years now (10+) with great
results. I figured it was pretty standard "least squares fit" type of
code. What I have discovered recently is that this code seems to produce
*much* better fits than Octave and Excel and I assume Matlab although I
haven't tested that. Anyway, I will try to put together some tests to
show that my port to Octave works as stated. I was actually hoping that
someone else would grab it and make it fit the coding conventions of Octave
as I am definitely an Octave beginner and not much of a Math expert
either. I feel confident that I can show that it is a much better
algorithm though. I will try to post some results this weekend.
Thanks,
-Todd
Post by John W. Eaton
| Hope this is the correct place...
|
| Please find a replacement for the current polyfit in Octave
attached. I
| believe it produces coefficients that are orders of magnitude better than
| the current Octave polyfit function. I was told to submit it to the
| mailing lists as a bug, hopefully this is the right place.
|
| Thanks
| -Todd
What is the origin of this code? Is it derived from some other
source? Also, is there some way to avoid the for loops and vectorize
it?
Before it could be included in Octave as a replacement for the current
polyfit function, it would need to better match the coding conventions
of the rest of Octave, have a doc string (could use the one from the
existing polyfit function), copyright statement, etc.
Also, I think we would want to some examples that show that it
produces correct results (tests would be nice) and some performance
numbers.
jwe
--
A VTOL Project in my spare time
http://www.hanfordsite.com/
(some images from WWJ code I have been working on)
http://picasaweb.google.com/tve9999/WorldWindJavaImages/photo#s5093635840816720994
Dmitri A. Sergatskov
2007-09-14 03:26:32 UTC
Permalink
I would be interested to see if this "improvement" still holds
against wpolyfit from octave-forge.
I suspect that all improvements comes from betterpolyfit
normalizing inputs.

There were few discussions on
the list of polyfit returning bad
results for input x going from 0 to 1e6 or so.

Perhaps polyfit should add the scaling the way wpolyfit.
Alternatively, in case that this improvement would violate
bug-for-bug compatibility with matlab, some wornings should
be added to the help file.

Sincerely,

Dmitri.
--
Post by Todd Elliott
Hi John,
I completely agree with your concerns. As I think I mentioned in my
original e-mail to you, I aquired this code in the Java language some
years ago. I believe it was "freeware" as in no-license-whatsoever from
the following (now defunct link)
<http://www.intergalact.com/java/polyfit/polyfit.html>
I have used the code in Java and C forms for many years now (10+) with great
results. I figured it was pretty standard "least squares fit" type of
code. What I have discovered recently is that this code seems to produce
*much* better fits than Octave and Excel and I assume Matlab although I
haven't tested that. Anyway, I will try to put together some tests to
show that my port to Octave works as stated. I was actually hoping that
someone else would grab it and make it fit the coding conventions of Octave
as I am definitely an Octave beginner and not much of a Math expert either.
I feel confident that I can show that it is a much better algorithm
though. I will try to post some results this weekend.
Thanks,
-Todd
Post by John W. Eaton
| Hope this is the correct place...
|
| Please find a replacement for the current polyfit in Octave attached.
I
Post by John W. Eaton
| believe it produces coefficients that are orders of magnitude better
than
Post by John W. Eaton
| the current Octave polyfit function. I was told to submit it to the
| mailing lists as a bug, hopefully this is the right place.
|
| Thanks
| -Todd
What is the origin of this code? Is it derived from some other
source? Also, is there some way to avoid the for loops and vectorize
it?
Before it could be included in Octave as a replacement for the current
polyfit function, it would need to better match the coding conventions
of the rest of Octave, have a doc string (could use the one from the
existing polyfit function), copyright statement, etc.
Also, I think we would want to some examples that show that it
produces correct results (tests would be nice) and some performance
numbers.
jwe
Todd Elliott
2007-09-14 22:04:44 UTC
Permalink
Hi Dmitri,
Thank you for your reply. I didn't even know about octave-forge. Now I
really feel out of my element :) That looks like a great resource. What
package is wpolyfit in? I didn't see it in the polynomial package. I
would like to include it in my testing. I'm starting to think that the
code I am wanting to submit might be better to include in something like
octave-forge to preserve compatibility with Matlab. I will say that I do
believe there may very well be a *bug* in the current Octave polyfit
function though. It might be something to do with a "sign" problem. I
noticed the coefficients it produces do not alternate like I would expect
(keep in mind that I do not know what I am doing). For example, when
trying to fit a thermistor response to temperature, I noticed that the
Octave polyfit function seems to produce much better results with a 4th
order fit vs 5th order fit. This is not the case with 'betterpolyfit' and
other thermistor fits I have observed. It seems as if there might be some
issue that arises with odd vs even order fits with the current polyfit
function.

I would like to show some results from the current polyfit, betterpolyfit,
and wpolyfit.
Please let me know what octave-forge package wpolyfit is in. My apologies
if you already mentioned it and I missed it.
Thanks,
-Todd
Post by Dmitri A. Sergatskov
I would be interested to see if this "improvement" still holds
against wpolyfit from octave-forge.
I suspect that all improvements comes from betterpolyfit
normalizing inputs.
There were few discussions on
the list of polyfit returning bad
results for input x going from 0 to 1e6 or so.
Perhaps polyfit should add the scaling the way wpolyfit.
Alternatively, in case that this improvement would violate
bug-for-bug compatibility with matlab, some wornings should
be added to the help file.
Sincerely,
Dmitri.
--
Post by Todd Elliott
Hi John,
I completely agree with your concerns. As I think I mentioned in my
original e-mail to you, I aquired this code in the Java language some
years ago. I believe it was "freeware" as in no-license-whatsoever
from
Post by Todd Elliott
the following (now defunct link)
<http://www.intergalact.com/java/polyfit/polyfit.html>
I have used the code in Java and C forms for many years now (10+) with
great
Post by Todd Elliott
results. I figured it was pretty standard "least squares fit" type of
code. What I have discovered recently is that this code seems to
produce
Post by Todd Elliott
*much* better fits than Octave and Excel and I assume Matlab although I
haven't tested that. Anyway, I will try to put together some tests
to
Post by Todd Elliott
show that my port to Octave works as stated. I was actually hoping
that
Post by Todd Elliott
someone else would grab it and make it fit the coding conventions of
Octave
Post by Todd Elliott
as I am definitely an Octave beginner and not much of a Math expert
either.
Post by Todd Elliott
I feel confident that I can show that it is a much better algorithm
though. I will try to post some results this weekend.
Thanks,
-Todd
Post by John W. Eaton
| Hope this is the correct place...
|
| Please find a replacement for the current polyfit in Octave
attached.
Post by Todd Elliott
I
Post by John W. Eaton
| believe it produces coefficients that are orders of magnitude better
than
Post by John W. Eaton
| the current Octave polyfit function. I was told to submit it to
the
Post by Todd Elliott
Post by John W. Eaton
| mailing lists as a bug, hopefully this is the right place.
|
| Thanks
| -Todd
What is the origin of this code? Is it derived from some other
source? Also, is there some way to avoid the for loops and vectorize
it?
Before it could be included in Octave as a replacement for the current
polyfit function, it would need to better match the coding conventions
of the rest of Octave, have a doc string (could use the one from the
existing polyfit function), copyright statement, etc.
Also, I think we would want to some examples that show that it
produces correct results (tests would be nice) and some performance
numbers.
jwe
--
A VTOL Project in my spare time
http://www.hanfordsite.com/
(some images from WWJ code I have been working on)
http://picasaweb.google.com/tve9999/WorldWindJavaImages/photo#s5093635840816720994
Dmitri A. Sergatskov
2007-09-15 01:32:07 UTC
Permalink
It is in "optim" (optimization) directory.
I missed what is your platform, but
octave-forge packages usually included with octave installs...

Dmitri.
Post by Todd Elliott
Hi Dmitri,
Thank you for your reply. I didn't even know about octave-forge. Now I
really feel out of my element :) That looks like a great resource. What
package is wpolyfit in? I didn't see it in the polynomial package. I
would like to include it in my testing. I'm starting to think that the
code I am wanting to submit might be better to include in something like
octave-forge to preserve compatibility with Matlab. I will say that I do
believe there may very well be a *bug* in the current Octave polyfit
function though. It might be something to do with a "sign" problem. I
noticed the coefficients it produces do not alternate like I would expect
(keep in mind that I do not know what I am doing). For example, when
trying to fit a thermistor response to temperature, I noticed that the
Octave polyfit function seems to produce much better results with a 4th
order fit vs 5th order fit. This is not the case with 'betterpolyfit' and
other thermistor fits I have observed. It seems as if there might be some
issue that arises with odd vs even order fits with the current polyfit
function.
I would like to show some results from the current polyfit, betterpolyfit,
and wpolyfit.
Please let me know what octave-forge package wpolyfit is in. My apologies
if you already mentioned it and I missed it.
Thanks,
-Todd
Post by Dmitri A. Sergatskov
I would be interested to see if this "improvement" still holds
against wpolyfit from octave-forge.
I suspect that all improvements comes from betterpolyfit
normalizing inputs.
There were few discussions on
the list of polyfit returning bad
results for input x going from 0 to 1e6 or so.
Perhaps polyfit should add the scaling the way wpolyfit.
Alternatively, in case that this improvement would violate
bug-for-bug compatibility with matlab, some wornings should
be added to the help file.
Sincerely,
Dmitri.
--
Post by Todd Elliott
Hi John,
I completely agree with your concerns. As I think I mentioned in my
original e-mail to you, I aquired this code in the Java language some
years ago. I believe it was "freeware" as in no-license-whatsoever
from
Post by Dmitri A. Sergatskov
Post by Todd Elliott
the following (now defunct link)
< http://www.intergalact.com/java/polyfit/polyfit.html>
I have used the code in Java and C forms for many years now (10+) with
great
Post by Dmitri A. Sergatskov
Post by Todd Elliott
results. I figured it was pretty standard "least squares fit" type of
code. What I have discovered recently is that this code seems to
produce
Post by Dmitri A. Sergatskov
Post by Todd Elliott
*much* better fits than Octave and Excel and I assume Matlab although I
haven't tested that. Anyway, I will try to put together some tests
to
Post by Dmitri A. Sergatskov
Post by Todd Elliott
show that my port to Octave works as stated. I was actually hoping
that
Post by Dmitri A. Sergatskov
Post by Todd Elliott
someone else would grab it and make it fit the coding conventions of
Octave
Post by Dmitri A. Sergatskov
Post by Todd Elliott
as I am definitely an Octave beginner and not much of a Math expert
either.
Post by Dmitri A. Sergatskov
Post by Todd Elliott
I feel confident that I can show that it is a much better algorithm
though. I will try to post some results this weekend.
Thanks,
-Todd
Post by John W. Eaton
| Hope this is the correct place...
|
| Please find a replacement for the current polyfit in Octave
attached.
Post by Dmitri A. Sergatskov
Post by Todd Elliott
I
Post by John W. Eaton
| believe it produces coefficients that are orders of magnitude better
than
Post by John W. Eaton
| the current Octave polyfit function. I was told to submit it to
the
Post by Dmitri A. Sergatskov
Post by Todd Elliott
Post by John W. Eaton
| mailing lists as a bug, hopefully this is the right place.
|
| Thanks
| -Todd
What is the origin of this code? Is it derived from some other
source? Also, is there some way to avoid the for loops and vectorize
it?
Before it could be included in Octave as a replacement for the current
polyfit function, it would need to better match the coding conventions
of the rest of Octave, have a doc string (could use the one from the
existing polyfit function), copyright statement, etc.
Also, I think we would want to some examples that show that it
produces correct results (tests would be nice) and some performance
numbers.
jwe
--
A VTOL Project in my spare time
http://www.hanfordsite.com/
(some images from WWJ code I have been working on)
http://picasaweb.google.com/tve9999/WorldWindJavaImages/photo#s5093635840816720994
Dmitri A. Sergatskov
2007-09-16 19:03:08 UTC
Permalink
To get decent results for polynomial fits of high order
one should really use set of orthogonal polynomial functions.
(Chebyshev's polynomials are popular choice.)
I assume that pivoting used by betterpolyfit somehow
helps with numerical instabilities, but it will break later.

Having said all this I would add that for this particular
data polynomial fit does not look like a good idea (the
data just does not look polynomial).
The "standard" way to describe thermistor output is
Steinhart-Hart equation (look it up) -- essentially
1/T as a polynomial of log(R).

Sincerely,

Dmitri.
Hi All,
Well, I had a chance to do a little bit of testing and the results were not
exactly what I expected, but betterpolyfit definitely came out on top so
far. There is definitely something wrong with the Octave polyfit and
Octaveforge wpolyfit functions. Before continuing, I thought I would
share some results. It seems that the coefficients produced by polyfit
and wpolyfit start breaking down for fits of 5th order and higher with
XRange values approching zero. Anybody have ideas?
Here is a pdf with some plots showing results
hanfordsite.com/testing/polyfit_testing.pdf
Regards,
-Todd
Post by Dmitri A. Sergatskov
It is in "optim" (optimization) directory.
I missed what is your platform, but
octave-forge packages usually included with octave installs...
Dmitri.
Dmitri A. Sergatskov
2007-09-16 19:54:41 UTC
Permalink
Post by Todd Elliott
Hi Dmitri,
Thanks for the info on the Steinhart equation. I have played around with
it, but I find that with betterpolyfit I can get a much better fit than with
the "standard" way of doing it. Using betterpolyfit, I can get the max
error down to 0.05 degrees. I have heard others say "that doesn't look
polynomial" or "too high of an order" before, but what does that mean?
If you can get a better fit over the entire range with a polynomial
equation, then why not?
Technically speaking you can plot n-1 degree polynomial
that will go exactly through n points, yet it would not
necessary make a good fit.

Try to plot (data - fit)
curve, or try to plot derivative of the fit -- I bet those
would look ugly (if you want to be more scientific --
plot a histogram of the residuals, -- does it look gaussian ?).
Post by Todd Elliott
Also, what do you mean it will break later? Can you give me an example?
No. What am saying is that fitting polynomial in the form of

sum by i of a_i * x^n-i

is numerically unstable for high n. While there are ways to
mitigate those instabilities, the better way is to fit
to a set of orthogonal functions.

But at the end it is all is what works for you. I just do not
think (at least not yet) that polyfit function in octave
have to be replaced.
Post by Todd Elliott
Thanks,
-Todd
Sincerely,

Dmitri.
--
Todd Elliott
2007-09-16 21:15:00 UTC
Permalink
Hi Dmitri and All,

I agree that polyfit may not need to be replaced, but from what I have seen
playing with it is that it probably a) needs to be fixed for fits of 5th
order and higher or b) limited to 4th order and lower
I believe there very well may be an actual bug in the implementation for 5th
order and higher fits.

The thermistor example is not all I have used the pivot-row based algorithm
for over the years. It always seems to do a good job for me. You are
right that for someone like me, the fact that it works is the most
important thing. btw, here is an example of a highly respected Scientific
measurement company that also used a 5th order polyfit for a thermistor
www.campbellsci.com/documents/manuals/107-lc.pdf They may have also used a
pivot-row based algorithm since their fit looks a lot like the fit
betterpolyfit would produce... or maybe polyfit and wpolyfit are broken for
a 5th order fit? In practice, I use a 5th order fit with a slightly
reduced temperature range compared to this to get a very nice fit with
minimal error due to the nature of a polynomial function. This can easily
be implemented in an 8-bit embedded system for highly accurate temperature
measurements using thermistors.

I have plotted the (data-fit) and the result as I mentioned in my original
post was that the betterpolyfit produced error several (3-4 times) orders of
magnitude less than the original polyfit function. As for plotting
residuals, and histograms, I'm afraid I will be leaving that for someone
else to do. I have an intuitive feeling that you will find the pivot-row
algorithm will come out on top every time.... possibly at the expense of a
few more computation cycles while generating your coefficients? (still
less than one second on my machine).

I have placed the temp_c, resistance, mv file here
http://www.hanfordsite.com/testing/thermistor.csv for reference in relation
to the polyfit_testing.pdf document.

Also, I have added some headers stating where the original source and
adaptation came from with the betterpolyfit function here
http://www.hanfordsite.com/testing/betterpolyfit.zip As stated in the
source, as far as I know and am concerned, there are no restrictions
whatsoever to this function.

Unless there is a compelling response to continue, I think I will probably
just leave this issue as it is. Hopefully someone else will also find
this function useful as I have.

Thanks,
Todd
Post by Dmitri A. Sergatskov
Post by Todd Elliott
Hi Dmitri,
Thanks for the info on the Steinhart equation. I have played around
with
Post by Todd Elliott
it, but I find that with betterpolyfit I can get a much better fit than
with
Post by Todd Elliott
the "standard" way of doing it. Using betterpolyfit, I can get the
max
Post by Todd Elliott
error down to 0.05 degrees. I have heard others say "that doesn't
look
Post by Todd Elliott
polynomial" or "too high of an order" before, but what does that mean?
If you can get a better fit over the entire range with a polynomial
equation, then why not?
Technically speaking you can plot n-1 degree polynomial
that will go exactly through n points, yet it would not
necessary make a good fit.
Try to plot (data - fit)
curve, or try to plot derivative of the fit -- I bet those
would look ugly (if you want to be more scientific --
plot a histogram of the residuals, -- does it look gaussian ?).
Post by Todd Elliott
Also, what do you mean it will break later? Can you give me an
example?
No. What am saying is that fitting polynomial in the form of
sum by i of a_i * x^n-i
is numerically unstable for high n. While there are ways to
mitigate those instabilities, the better way is to fit
to a set of orthogonal functions.
But at the end it is all is what works for you. I just do not
think (at least not yet) that polyfit function in octave
have to be replaced.
Post by Todd Elliott
Thanks,
-Todd
Sincerely,
Dmitri.
--
--
A VTOL Project in my spare time
http://www.hanfordsite.com/
(some images from WWJ code I have been working on)
http://picasaweb.google.com/tve9999/WorldWindJavaImages/photo#s5093635840816720994
Todd Elliott
2007-09-16 21:21:13 UTC
Permalink
Just to be clear, I meant 3-4 orders of magnitude, not 3-4 times. i.e.
a factor of 10,000
Post by Todd Elliott
I have plotted the (data-fit) and the result as I mentioned in my original
post was that the betterpolyfit >>produced error several (3-4 times) orders
of magnitude less than the original polyfit function. As for plotting

-Todd

Todd Elliott
2007-09-16 19:17:09 UTC
Permalink
Hi Dmitri,

Thanks for the info on the Steinhart equation. I have played around with
it, but I find that with betterpolyfit I can get a much better fit than with
the "standard" way of doing it. Using betterpolyfit, I can get the max
error down to 0.05 degrees. I have heard others say "that doesn't look
polynomial" or "too high of an order" before, but what does that mean?
If you can get a better fit over the entire range with a polynomial
equation, then why not?

Also, what do you mean it will break later? Can you give me an example?

Thanks,
-Todd
Post by Dmitri A. Sergatskov
To get decent results for polynomial fits of high order
one should really use set of orthogonal polynomial functions.
(Chebyshev's polynomials are popular choice.)
I assume that pivoting used by betterpolyfit somehow
helps with numerical instabilities, but it will break later.
Having said all this I would add that for this particular
data polynomial fit does not look like a good idea (the
data just does not look polynomial).
The "standard" way to describe thermistor output is
Steinhart-Hart equation (look it up) -- essentially
1/T as a polynomial of log(R).
Sincerely,
Dmitri.
Hi All,
Well, I had a chance to do a little bit of testing and the results were
not
exactly what I expected, but betterpolyfit definitely came out on top
so
far. There is definitely something wrong with the Octave polyfit and
Octaveforge wpolyfit functions. Before continuing, I thought I
would
share some results. It seems that the coefficients produced by
polyfit
and wpolyfit start breaking down for fits of 5th order and higher with
XRange values approching zero. Anybody have ideas?
Here is a pdf with some plots showing results
hanfordsite.com/testing/polyfit_testing.pdf
Regards,
-Todd
Post by Dmitri A. Sergatskov
It is in "optim" (optimization) directory.
I missed what is your platform, but
octave-forge packages usually included with octave installs...
Dmitri.
--
A VTOL Project in my spare time
http://www.hanfordsite.com/
(some images from WWJ code I have been working on)
http://picasaweb.google.com/tve9999/WorldWindJavaImages/photo#s5093635840816720994
Todd Elliott
2007-09-16 18:09:31 UTC
Permalink
Hi All,

Well, I had a chance to do a little bit of testing and the results were not
exactly what I expected, but betterpolyfit definitely came out on top so
far. There is definitely something wrong with the Octave polyfit and
Octaveforge wpolyfit functions. Before continuing, I thought I would
share some results. It seems that the coefficients produced by polyfit
and wpolyfit start breaking down for fits of 5th order and higher with
XRange values approching zero. Anybody have ideas?

Here is a pdf with some plots showing results
hanfordsite.com/testing/polyfit_testing.pdf

Regards,
-Todd
Post by Dmitri A. Sergatskov
It is in "optim" (optimization) directory.
I missed what is your platform, but
octave-forge packages usually included with octave installs...
Dmitri.
Loading...