Solved

# Interpolating unknows from polynomial curve

Posted on 2014-12-17
155 Views
Hi,

I'm am trying to script an existing workflow and I am struggling with the last part.

The users have a standard curve of points as part of the data (x = percentage, y=value) and in an external software (Prism, GraphPad) utilise a 6th order polynomial non-linear regression to utilise that curve and then interpolate unknowns from that.

I have my x and y in 2 arrays ready to go.
In trying to recreate this I have been having a try with Math::Polynomial but the interpolate function of that seems to be using the curve data to predict the polynomial needed so I am a little confused.

Does anyone have a hint as to how best to approach this in perl by reading in curve data then reading out the unknowns?

Thanks
0
Question by:Epipgx
• 3
• 2

LVL 84

Expert Comment

If you give Math::Polynomial->interpolate 7 data points, it can generate a 6th order polynomial that passes through those points, you can then evaluate that polynomial at other points.
If that's not what you want, perhaps you want to use Algorithm::CurveFit or PDL::Fit::Polynomial
0

Author Comment

Thanks Ozo,

So there isn't a way to specify the the nth order? It is (as reading a bit more on the subject I think it should be) driven from the curve itself and produces the best order fit for your data?
If I have understood that correctly can you point me to an example where once the curve is fit using that value to read unknowns from it?
0

LVL 84

Accepted Solution

ozo earned 500 total points
#!/usr/bin/perl
package  Epipgx;
use strict;
use warnings;
use Carp;
use Math::MatrixReal;
sub new{
my(\$self,\$o) = @_;
my \$class = ref \$self || \$self;
\$o//=6;
croak "order must be positive" if \$o<0;
return bless {order=>\$o},\$class;
}
sub interpolate {
my(\$self, \$xvalues, \$yvalues) = @_;
my \$o=\$self->{order}//6;
\$o=\$#\$xvalues if \$o>\$#\$xvalues;
my @ix=map{my\$x=\$_;[map{\$x**\$_}0..\$o]}@\$xvalues;
my @iy=@\$yvalues;
my @XtX;
for my \$x( @ix ){
for my \$i ( 0..\$o ){
for my \$j ( 0..\$o ){
\$XtX[\$i][\$j] += \$x->[\$i]*\$x->[\$j];
}
}
}
my \$XtX = Math::MatrixReal->new_from_rows( [@XtX] );
my \$inv = \$XtX->inverse() or warn \$!;
my \$iY = Math::MatrixReal->new_from_cols( [[@iy]] );
my \$iX = Math::MatrixReal->new_from_rows( [@ix] );
my \$P = \$inv * ~\$iX;
my \$W = \$P * \$iY;
\$self->{a}=[map{\$W->element(\$_+1,1)}0..\$o];
}
sub evaluate{
my (\$self, \$x) = @_;
my \$y=0;
my \$X=1;
for my \$w( @{\$self->{a}} ){
\$y+=\$w*\$X;
\$X*=\$x;
}
return \$y;
}
sub _stringify{
my(\$self)=@_;
return join(" + ",map{(\$self->{a}[\$_]//0).(\$_>0&&"*x").(\$_>1&&"^\$_")}reverse 0..\$self->{order})=~s/\+ -/- /gr;
}
1;

Epipgx - Perl class for Interpolating unknows from polynomial curve

use Epipgx;

my \$p=Epipgx->new(6);   # 6th order polynomial

\$p->interpolate([0..8], [0, 1, 1, 2, 3, 5, 8, 13, 21]);

print "\$p\n";

\$y = \$p->evaluate(0.5);

=cut
0

Author Comment

Thanks Ozo that is fantastic, you have gone way beyond for me there and is much appreciated.
0

Author Closing Comment

Fantastically thorough answer, I was hoping for guidance but received a fully formed solution to my question, brilliant
0

## Featured Post

### Suggested Solutions

Relative Frequency Distribution 4 24
Probability Calculation 2 42
Independent Events 4 50
Exam question 48 99
On Microsoft Windows, if  when you click or type the name of a .pl file, you get an error "is not recognized as an internal or external command, operable program or batch file", then this means you do not have the .pl file extension associated with …
Lithium-ion batteries area cornerstone of today's portable electronic devices, and even though they are relied upon heavily, their chemistry and origin are not of common knowledge. This article is about a device on which every smartphone, laptop, an…
Explain concepts important to validation of email addresses with regular expressions. Applies to most languages/tools that uses regular expressions. Consider email address RFCs: Look at HTML5 form input element (with type=email) regex pattern: T…
This video gives you a great overview about bandwidth monitoring with SNMP and WMI with our network monitoring solution PRTG Network Monitor (https://www.paessler.com/prtg). If you're looking for how to monitor bandwidth using netflow or packet s…