polyfit.php 1.8 KB

12345678910111213141516171819202122232425262728293031323334353637383940414243444546474849505152535455565758596061626364656667686970717273
  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. function printpoly( $C = null ) {
  30. for($i = $C->m - 1; $i >= 0; --$i) {
  31. $r = $C->get($i, 0);
  32. if ( abs($r) <= pow(10, -9) )
  33. $r = 0;
  34. if ($i == $C->m - 1)
  35. echo $r . "x<sup>$i</sup>";
  36. else if ($i < $C->m - 1)
  37. echo " + " . $r . "x<sup>$i</sup>";
  38. else if ($i == 0)
  39. echo " + " . $r;
  40. }
  41. }
  42. $X = array(0,1,2,3,4,5);
  43. $Y = array(4,3,12,67,228, 579);
  44. $points = new Matrix(array($X, $Y));
  45. $points->toHTML();
  46. printpoly(polyfit($X, $Y, 4));
  47. echo '<hr />';
  48. $X = array(0,1,2,3,4,5);
  49. $Y = array(1,2,5,10,17, 26);
  50. $points = new Matrix(array($X, $Y));
  51. $points->toHTML();
  52. printpoly(polyfit($X, $Y, 2));
  53. echo '<hr />';
  54. $X = array(0,1,2,3,4,5,6);
  55. $Y = array(-90,-104,-178,-252,-26, 1160, 4446);
  56. $points = new Matrix(array($X, $Y));
  57. $points->toHTML();
  58. printpoly(polyfit($X, $Y, 5));
  59. echo '<hr />';
  60. $X = array(0,1,2,3,4);
  61. $Y = array(mt_rand(0, 10), mt_rand(40, 80), mt_rand(240, 400), mt_rand(1800, 2215), mt_rand(8000, 9000));
  62. $points = new Matrix(array($X, $Y));
  63. $points->toHTML();
  64. printpoly(polyfit($X, $Y, 3));
  65. ?>