ILIAS  Release_4_0_x_branch Revision 61816
 All Data Structures Namespaces Files Functions Variables Groups Pages
polyfit.php
Go to the documentation of this file.
1 <?php
2 require_once "../Matrix.php";
3 /*
4 * @package JAMA
5 * @author Michael Bommarito
6 * @author Paul Meagher
7 * @version 0.1
8 *
9 * Function to fit an order n polynomial function through
10 * a series of x-y data points using least squares.
11 *
12 * @param $X array x values
13 * @param $Y array y values
14 * @param $n int order of polynomial to be used for fitting
15 * @returns array $coeffs of polynomial coefficients
16 * Pre-Conditions: the system is not underdetermined: sizeof($X) > $n+1
17 */
18 function polyfit($X, $Y, $n) {
19  for ($i = 0; $i < sizeof($X); ++$i)
20  for ($j = 0; $j <= $n; ++$j)
21  $A[$i][$j] = pow($X[$i], $j);
22  for ($i=0; $i < sizeof($Y); ++$i)
23  $B[$i] = array($Y[$i]);
24  $matrixA = new Matrix($A);
25  $matrixB = new Matrix($B);
26  $C = $matrixA->solve($matrixB);
27  return $C->getMatrix(0, $n, 0, 1);
28 }
29 
30 function printpoly( $C = null ) {
31  for($i = $C->m - 1; $i >= 0; --$i) {
32  $r = $C->get($i, 0);
33  if ( abs($r) <= pow(10, -9) )
34  $r = 0;
35  if ($i == $C->m - 1)
36  echo $r . "x<sup>$i</sup>";
37  else if ($i < $C->m - 1)
38  echo " + " . $r . "x<sup>$i</sup>";
39  else if ($i == 0)
40  echo " + " . $r;
41  }
42 }
43 
44 $X = array(0,1,2,3,4,5);
45 $Y = array(4,3,12,67,228, 579);
46 $points = new Matrix(array($X, $Y));
47 $points->toHTML();
48 printpoly(polyfit($X, $Y, 4));
49 
50 echo '<hr />';
51 
52 $X = array(0,1,2,3,4,5);
53 $Y = array(1,2,5,10,17, 26);
54 $points = new Matrix(array($X, $Y));
55 $points->toHTML();
56 printpoly(polyfit($X, $Y, 2));
57 
58 echo '<hr />';
59 
60 $X = array(0,1,2,3,4,5,6);
61 $Y = array(-90,-104,-178,-252,-26, 1160, 4446);
62 $points = new Matrix(array($X, $Y));
63 $points->toHTML();
64 printpoly(polyfit($X, $Y, 5));
65 
66 echo '<hr />';
67 
68 $X = array(0,1,2,3,4);
69 $Y = array(mt_rand(0, 10), mt_rand(40, 80), mt_rand(240, 400), mt_rand(1800, 2215), mt_rand(8000, 9000));
70 $points = new Matrix(array($X, $Y));
71 $points->toHTML();
72 printpoly(polyfit($X, $Y, 3));
73 ?>