* $data = array(array($x1, $y1), array($x2, $y2));
* $linear = new zc_math_linear($data);
* $linear->calculate();
* // y = k*x + m
* echo $linear->k;
* echo $linear->m;
*
*
* Linear regression formula found at {@link http://en.wikipedia.org/wiki/Numerical_methods_for_linear_least_squares#Computation }
*
* @author Han Lin Yap < http://zencodez.net/ >
* @copyright 2010 zencodez.net
* @license http://creativecommons.org/licenses/by-sa/3.0/
* @package math
* @version 1.0 - 2010-07-29
*/
class zc_math_linear {
/**
* Contains the coordinates
* @var array
*/
public $data;
// needed for linear regression formula
public $sum_x;
public $sum_x2;
public $sum_y;
public $sum_xy;
public $number_of_coords;
// y = k*x + m
/**
* @var float
*/
public $k;
/**
* @var float
*/
public $m;
/**
* @param array $data
*/
function __construct($data) {
$this->data = $data;
}
/**
* Calculate and make formula
* To understand better, {@link http://en.wikipedia.org/wiki/Numerical_methods_for_linear_least_squares#Computation see formula}
* @since 1.0
*/
public function calculate() {
$this->number_of_coords = count($this->data);
$this->sum_x = $this->sum_by_formula($this->data, '$x[0]');
$this->sum_x2 = $this->sum_by_formula($this->data, 'pow($x[0],2)');
$this->sum_y = $this->sum_by_formula($this->data, '$x[1]');
$this->sum_xy = $this->sum_by_formula($this->data, '$x[0] * $x[1]');
$this->k = ($this->number_of_coords * $this->sum_xy - $this->sum_y * $this->sum_x) / ($this->number_of_coords * $this->sum_x2 - pow($this->sum_x, 2));
$this->m = ($this->sum_y * $this->sum_x2 - $this->sum_x * $this->sum_xy) / ($this->number_of_coords * $this->sum_x2 - pow($this->sum_x, 2));
}
/**
* Gets all by a formula
* Turn a multidimention array to a ordinary array with values
* @param array $data Multidimention array
* @param string $formula How to get what
* @return array Get values in multidimention by formula
* @since 1.0
*/
protected function _get_values_by_formula($data, $formula) {
return array_map(create_function('$x', 'return ' . $formula . ';'), $data);
}
/**
* Sum all by a formula
* @param array $data Data to sum
* @param string $formula Formula to sum
* @return float sum of $data by formula
* @since 1.0
*/
public function sum_by_formula($data, $formula) {
return array_sum($this->_get_values_by_formula($data, $formula));
}
}
/**
* Polynomial regression
*
* Example usage
*
* $polynomial = new zc_math_polynomial($data);
* echo print_r($polynomial->calculate_coefficients());
* echo "![](" . $polynomial->render_formula() . ")
";
* echo "" . $polynomial->render_formula('') . "";
* echo $polynomial->get_y(1200);
*
* {@link http://www.trentfguidry.net/post/2009/06/30/Matrix-Crout-LU-decomposition.aspx }
* {@link http://www.trentfguidry.net/post/2009/07/19/Linear-multiple-regression.aspx }
* {@link http://www.trentfguidry.net/post/2009/08/01/Linear-regression-polynomial-coefficients.aspx }
* @author Han Lin Yap < http://zencodez.net/ >
* @copyright 2010 zencodez.net
* @license http://creativecommons.org/licenses/by-sa/3.0/
* @package math
* @version 1.0 - 2010-07-30
*/
class zc_math_polynomial extends zc_math_linear {
/**
* The calculated formula in array
* @var array
*/
public $result;
/**
* Order to calculate
* @var integer
*/
public $order;
/**
* Initialize
* @param array $data Format array(array($x1, $y2), array($x2, $y2))
* @param integer $order optional What order the formula should be
* @since 1.0
*/
function __construct($data, $order=4) {
$this->data = $data;
$this->set_order($order);
}
/**
* Change order
* @param integer $order What order the formula should be
* @since 1.0
*/
public function set_order($order) {
$this->order = $order;
}
/**
* Find the formula
* @param integer $order optional What order the formula should be
* @return array The found formula in array
* @since 1.0
*/
public function calculate_coefficients($order=false) {
if ($order)
$this->set_order($order);
$x = $this->_get_values_by_formula($this->data, '$x[0]');
$y = $this->_get_values_by_formula($this->data, '$x[1]');
$z = array_fill(0, count($y), array());
for ($i = 0 ; $i < count($y) ; $i++) {
for ($j = 0 ; $j <= $this->order ; $j++) {
$z[$i][$j] = pow($x[$i],$j);
}
}
$this->result = $this->_polynomial_regress($z, $y);
return $this->result;
}
/**
* Render formula for wolframalpha, google chart, url or plain math text
*
* NOTE! {@link zc_math_polynomial::calculate_coefficients()} before calling this method
* Example usage
*
* render_formula()
* render_formula('wolframalpha')
* render_formula('')
* render_formula('text');
*
*
* @param string $type optional Available types wolframalpha, text, google or ''
* @param string|boolean $replace_x optional If you want to customize x output
* @return string URL or text, depends on $type
* @uses zc_math_polynomial::calculate_coefficients() Call it before to have something to render
* @since 1.0
*/
public function render_formula($type='google', $replace_x=false) {
$formula = array();
$text = '';
if ($type=='text')
$text = '*';
if (!$replace_x)
$x = 'x';
else
$x = $replace_x;
foreach($this->result AS $k => $v) {
if ($k == 0)
$formula[] = $v[0];
elseif ($k == 1)
$formula[] = $v[0] . $text . $x;
elseif ($type!='google')
$formula[] = $v[0] . $text . $x . "^". $k;
else
$formula[] = $v[0] . "x^{". $k . "}";
}
if ($type=='wolframalpha')
return 'http://www.wolframalpha.com/input/?i=' . implode("+%2b+", $formula);
elseif ($type=='text')
return implode("%2b", $formula);
elseif ($type!='google')
return 'y = ' . implode(" + ", $formula);
else
return 'http://chart.apis.google.com/chart?cht=tx&chl=y+=+' . implode("%2b", $formula);
}
/**
* Get y value from x coordinate
*
* NOTE! {@link zc_math_polynomial::calculate_coefficients()} before calling this method
*
* @param integer $x X-coordinate
* @return integer Y-coordinate
* @uses zc_math_polynomial::calculate_coefficients() Call it before to have something to check on
* @since 1.0
*/
public function get_y($x) {
$sum = 0;
foreach($this->result AS $k => $v) {
if ($k == 0)
$sum += $v[0];
else
$sum += $v[0] * pow($x,$k);
}
return $sum;
}
private function _polynomial_regress($z, $y) {
$z_transpose = zc_math::matrix_transpose($z);
$l = zc_math::matrix_multiply($z_transpose, $z);
$r = zc_math::matrix_multiply($z_transpose, zc_math::matrix_transpose(array($y)));
return $this->_solve_for($l, $r);
}
private function _solve_for($l, $r) {
$resultMatrix = array_fill(0, count($l[0]), array());
$resDecomp = $this->_lu_decompose($l);
$nP = $resDecomp['PivotArray'];
$lMatrix = $resDecomp['L'];
$uMatrix = $resDecomp['U'];
for ($k = 0 ; $k < count($r[0]) ; $k++) {
$sum = 0.0;
$dMatrix = array_fill(0, count($l), array());
$dMatrix[0][0] = $r[$nP[0]][$k] / $lMatrix[0][0];
for ($i = 1; $i < count($l); $i++) {
$sum = 0.0;
for ($j = 0; $j < $i; $j++) {
$sum += $lMatrix[$i][$j] * $dMatrix[$j][0];
}
$dMatrix[$i][0] = ($r[$nP[$i]][$k] - $sum) / $lMatrix[$i][$i];
}
$resultMatrix[count($l) - 1][$k] = $dMatrix[count($l) - 1][0];
for ($i = count($l) - 2; $i >= 0; $i--) {
$sum = 0.0;
for ($j = $i + 1; $j < count($l); $j++) {
$sum += $uMatrix[$i][$j] * $resultMatrix[$j][$k];
}
$resultMatrix[$i][$k] = $dMatrix[$i][0] - $sum;
}
}
return $resultMatrix;
}
private function _lu_decompose($l) {
if (!defined('PHP_INT_MIN')) {
define('PHP_INT_MIN', ~PHP_INT_MAX);
}
$_rowCount = count($l);
$_columnCount = count($l[0]);
$uMatrix = array_fill(0, $_rowCount, array());
$lMatrix = array_fill(0, $_rowCount, array());
$workingUMatrix = $l;
$workingLMatrix = array_fill(0, $_rowCount, array());
$pivotArray = range(0, $_rowCount-1);
for ($i = 0; $i < $_rowCount; $i++) {
$maxRowRatio = PHP_INT_MIN;
$maxRow = -1;
$maxPosition = -1;
for ($j = $i; $j < $_rowCount; $j++) {
$rowSum = 0.0;
for ($k = $i; $k < $_columnCount; $k++) {
$rowSum += abs($workingUMatrix[$pivotArray[$j]][$k]);
}
$dCurrentRatio = abs($workingUMatrix[$pivotArray[$j]][$i]) / $rowSum;
if ($dCurrentRatio > $maxRowRatio) {
$maxRowRatio = abs($workingUMatrix[$pivotArray[$j]][$i] / $rowSum);
$maxRow = $pivotArray[$j];
$maxPosition = $j;
}
}
if ($maxRow != $pivotArray[$i]) {
$hold = $pivotArray[$i];
$pivotArray[$i] = $maxRow;
$pivotArray[$maxPosition] = $hold;
}
$rowFirstElementValue = $workingUMatrix[$pivotArray[$i]][$i];
for ($j = 0; $j < $_columnCount; $j++) {
if ($j < $i) {
$workingUMatrix[$pivotArray[$i]][$j] = 0.0;
} elseif ($j == $i) {
$workingLMatrix[$pivotArray[$i]][$j] = $rowFirstElementValue;
$workingUMatrix[$pivotArray[$i]][$j] = 1.0;
} else {
$workingUMatrix[$pivotArray[$i]][$j] /= $rowFirstElementValue;
$workingLMatrix[$pivotArray[$i]][$j] = 0.0;
}
}
for ($k = $i + 1; $k < $_rowCount; $k++) {
$rowFirstElementValue = $workingUMatrix[$pivotArray[$k]][$i];
for ($j = 0; $j < $_rowCount; $j++) {
if ($j < $i) {
$workingUMatrix[$pivotArray[$k]][$j] = 0.0;
} elseif ($j == $i) {
$workingLMatrix[$pivotArray[$k]][$j] = $rowFirstElementValue;
$workingUMatrix[$pivotArray[$k]][$j] = 0.0;
} else {
$workingUMatrix[$pivotArray[$k]][$j] = $workingUMatrix[$pivotArray[$k]][$j] - $rowFirstElementValue * $workingUMatrix[$pivotArray[$i]][$j];
}
}
}
}
for ($i = 0; $i < $_rowCount; $i++) {
for ($j = 0; $j < $_rowCount; $j++) {
$uMatrix[$i][$j] = $workingUMatrix[$pivotArray[$i]][$j];
$lMatrix[$i][$j] = $workingLMatrix[$pivotArray[$i]][$j];
}
}
return array('U' => $uMatrix, 'L' => $lMatrix, 'PivotArray' => $pivotArray);
}
}
?>