Polynomial regression: Difference between revisions

This script depends on the <tt>Math::MatrixReal</tt> CPAN module to compute matrix determinants.

<lang Perl>

## CAUTION - THE BELOW SCRIPT DOES NOT WORK. Example:-

# $ perl polyfit.pl

# Please enter the values for the x coordinates, each delimited by a space. (Ex: 0 1 2 3)

# 0 1 2 3 4 5 6 7 8 9

# Please enter the values for the y coordinates, each delimited by a space. (Ex: 0 1 2 3)

# 2.7 2.8 31.4 38.1 58.0 76.2 100.5 130.0 149.3 180.0

# An approximating polynomial for your dataset is -1x^7 + 12x^6 + -71x^5 + 259x^4 + -555x^3 + 630x^2 + -274x + 3.

##

# perl -e 'foreach $x (0..9) {print -1*$x**7 + 12*$x**6 + -71*$x**5 + 259*$x**4 + -555*$x**3 + 630*$x**2 + -274*$x + 3; print " ";}'

# 3 3 47 153 -165 -5617 -35337 -144603 -463117 -1254885

## (above is obviously wrong) - below is what it should be:

# perl -e 'foreach $x (0..9) {print 1.0848*$x**2+10.3552*$x-0.6164; print " ";}'

# -0.6164 10.8236 24.4332 40.2124 58.1612 78.2796 100.5676 125.0252 151.6524 180.4492

## (working polynomial above is from another implementation on this page)

##

use strict;

use feature 'say';

#This is a script to calculate an equation for a given set of coordinates.

#Input will be taken in sets of x and y. It can handle a grand total of 26 pairs.

#For matrix functions, we depend on the Math::MatrixReal package.

use Math::MatrixReal;

=pod

Step 1: Get each x coordinate all at once (delimited by " ") and each for y at once

on the next prompt in the same format (delimited by " ").

=cut

sub getPairs() {

my $buffer = <STDIN>;

chomp($buffer);

return split(" ", $buffer);

}

say("Please enter the values for the x coordinates, each delimited by a space. \(Ex: 0 1 2 3\)");

my @x = getPairs();

say("Please enter the values for the y coordinates, each delimited by a space. \(Ex: 0 1 2 3\)");

my @y = getPairs();

#This whole thing depends on the number of x's being the same as the number of y's

my $pairs = scalar(@x);

=pod

Step 2: Devise the base equation of our polynomial using the following idea

There is some polynomial of degree n (n == number of pairs - 1) such that

f(x)=ax^n + bx^(n-1) + ... yx + z

=cut

#Create an array of coefficients and their degrees with the format ("coefficent degree")

my @alphabet;

my @degrees;

for(my $alpha = "a", my $degree = $pairs - 1; $degree >= 0; $degree--, $alpha++) {

push(@alphabet, "$alpha");

push(@degrees, "$degree");

}

=pod

Step 3: Using the array of coeffs and their degrees, set up individual equations solving for

each coordinate pair. Why put it in this format? It interfaces witht he Math::MatrixReal package better this way.

=cut

my @coeffs;

for(my $count = 0; $count < $pairs; $count++) {

my $buffer = "[ ";

foreach (@degrees) {

$buffer .= (($x[$count] ** $_) . " ");

}

push(@coeffs, ($buffer . "]"));

}

my $row;

foreach (@coeffs) {

$row .= ("$_\n");

}

=pod

Step 4: We now have rows of x's raised to powers. With this in mind, we create a coefficient matrix.

=cut

my $matrix = Math::MatrixReal->new_from_string($row);

my $buffMatrix = $matrix->new_from_string($row);

=pod

Step 5: Now that we've gotten the matrix to do what we want it to do, we need to calculate the various determinants of the matrices

=cut

my $coeffDet = $matrix->det();

=pod

Step 6: Now that we have the determinant of the coefficient matrix, we need to find the determinants of the coefficient matrix with each column (1 at a time) replaced with the y values.

=cut

#NOTE: Unlike in Perl, matrix indices start at 1, not 0.

for(my $rows = my $column = 1; $column <= $pairs; $column++) {

#Reassign the values in the current column to the y values

foreach (@y) {

$buffMatrix->assign($rows, $column, $_);

$rows++;

}

#Find the values for the variables a, b, ... y, z in the original polynomial

#To round the difference of the determinants, I had to get creative

my $buffDet = $buffMatrix->det() / $coeffDet;

my $tempDet = int(abs($buffDet) + .5);

$alphabet[$column - 1] = $buffDet >= 0 ? $tempDet : 0 - $tempDet;

#Reset the buffer matrix and the row counter

$buffMatrix = $matrix->new_from_string($row);

$rows = 1;

}

=pod

Step 7: Now that we've found the values of a, b, ... y, z of the original polynomial, it's time to form our polynomial!

=cut

my $polynomial;

for(my $i = 0; $i < $pairs-1; $i++) {

if($alphabet[$i] == 0) {

next;

}

if($alphabet[$i] == 1) {

$polynomial .= ($degrees[$i] . " + ");

}

if($degrees[$i] == 1) {

$polynomial .= ($alphabet[$i] . "x" . " + ");

}

else {

$polynomial .= ($alphabet[$i] . "x^" . $degrees[$i] . " + ");

}

}

#Now for the last piece of the poly: the y-intercept.

$polynomial .= $alphabet[scalar(@alphabet)-1];

print("An approximating polynomial for your dataset is $polynomial.\n");

</lang>

<pre>Please enter the values for the x coordinates, each delimited by a space. (Ex: 0 1 2 3)

0 1 2 3 4 5 6 7 8 9 10

Please enter the values for the y coordinates, each delimited by a space. (Ex: 0 1 2 3)

1 6 17 34 57 86 121 162 209 262 321

An approximating polynomial for your dataset is 3x^2 + 2x + 1.</pre>