出典：http://PolynomialRegression.drque.net/
最新版は上記から取得してください。
PolynomialRegression.php のソース
class PolynomialRegression
{
private $xPowers;
private $xyPowers;
private $numberOfCoefficient;
private $forcedValue;
public function __construct( $numberOfCoefficient = 3 )
{
$this->numberOfCoefficient = $numberOfCoefficient;
$this->reset();
} // __construct
public function reset()
{
$this->forcedValue = array();
$this->xPowers = array();
$this->xyPowers = array();
$squares = ( $this->numberOfCoefficient - 1 ) * 2;
// Initialize power arrays.
for ( $index = 0; $index <= $squares; ++$index )
{
$this->xPowers[ $index ] = 0;
$this->xyPowers[ $index ] = 0;
}
} // reset
public function setDegree( $numberOfCoefficient )
{
$this->numberOfCoefficient = $numberOfCoefficient;
} // setDegree
public function setNumberOfCoefficient( $numberOfCoefficient )
{
$this->numberOfCoefficient = $numberOfCoefficient;
} // setNumberOfCoefficient
public function getNumberOfCoefficient( $numberOfCoefficient )
{
return $this->numberOfCoefficient;
} // getnumberOfCoefficient
public function setForcedCoefficient( $coefficient, $value )
{
$this->forcedValue[ $coefficient ] = $value;
} // setForcedCoefficient
public function getForcedCoefficient( $coefficient, $value )
{
$result = null;
if ( isset( $this->forcedValue[ $coefficient ] ) )
$result = $this->forcedValue[ $coefficient ];
return $result;
} // getForcedCoefficient
public function addData( $x, $y )
{
$squares = ( $this->numberOfCoefficient - 1 ) * 2;
// Remove the effect of the forced coefficient from this value.
foreach ( $this->forcedValue as $coefficient => $value )
{
$sub = bcpow( $x, $coefficient );
$sub = bcmul( $sub, $value );
$y = bcsub( $y, $sub );
}
// Accumulate new data to power sums.
for ( $index = 0; $index <= $squares; ++$index )
{
$this->xPowers[ $index ] =
bcadd( $this->xPowers[ $index ], bcpow( $x, $index ) );
$this->xyPowers[ $index ] =
bcadd
(
$this->xyPowers[ $index ],
bcmul( $y, bcpow( $x, $index ) )
);
}
} // addData
public function getCoefficients( $numberOfCoefficient = -1 )
{
// If no number of coefficients specified, use standard.
if ( $numberOfCoefficient == -1 )
$numberOfCoefficient = $this->numberOfCoefficient;
$matrix = array();
for ( $row = 0; $row < $numberOfCoefficient; ++$row )
{
$matrix[ $row ] = array();
for ( $column = 0; $column < $numberOfCoefficient; ++$column )
$matrix[ $row ][ $column ] =
$this->xPowers[ $row + $column ];
}
// Create augmented matrix by adding X*Y powers.
for ( $row = 0; $row < $numberOfCoefficient; ++$row )
$matrix[ $row ][ $numberOfCoefficient ] = $this->xyPowers[ $row ];
foreach ( $this->forcedValue as $coefficient => $value )
{
for ( $index = 0; $index < $numberOfCoefficient; ++$index )
{
$matrix[ $index ][ $coefficient ] = "0";
$matrix[ $coefficient ][ $index ] = "0";
}
$matrix[ $coefficient ][ $coefficient ] = "1";
$matrix[ $coefficient ][ $numberOfCoefficient ] = $value;
}
// Determine number of rows in matrix.
$rows = count( $matrix );
// Initialize done.
$isDone = array();
for ( $column = 0; $column < $rows; ++$column )
$isDone[ $column ] = false;
$order = array();
for ( $column = 0; $column < $rows; ++$column )
{
// Find a row to work with.
// A row that has a term in this column, and has not yet been
// reduced.
$activeRow = 0;
while ( ( ( 0 == $matrix[ $activeRow ][ $column ] )
|| ( $isDone[ $activeRow ] ) )
&& ( $activeRow < $rows ) )
{
++$activeRow;
}
// Do we have a term in this row?
if ( $activeRow < $rows )
{
// Remember the order.
$order[ $column ] = $activeRow;
// Normalize row--results in the first term being 1.
$firstTerm = $matrix[ $activeRow ][ $column ];
for ( $subColumn = $column; $subColumn <= $rows; ++$subColumn )
$matrix[ $activeRow ][ $subColumn ] =
bcdiv( $matrix[ $activeRow ][ $subColumn ], $firstTerm );
// This row is finished.
$isDone[ $activeRow ] = true;
// Subtract the active row from all rows that are not finished.
for ( $row = 0; $row < $rows; ++$row )
if ( ( ! $isDone[ $row ] )
&& ( 0 != $matrix[ $row ][ $column ] ) )
{
// Get first term in row.
$firstTerm = $matrix[ $row ][ $column ];
for ( $subColumn = $column; $subColumn <= $rows; ++$subColumn )
{
$accumulator = bcmul( $firstTerm, $matrix[ $activeRow ][ $subColumn ] );
$matrix[ $row ][ $subColumn ] =
bcsub( $matrix[ $row ][ $subColumn ], $accumulator );
}
}
}
}
// Reset done.
for ( $row = 0; $row < $rows; ++$row )
$isDone[ $row ] = false;
$coefficients = array();
for ( $column = ( $rows - 1 ); $column >= 0; --$column )
{
// The active row is based on order.
$activeRow = $order[ $column ];
// The active row is now finished.
$isDone[ $activeRow ] = true;
// For all rows not finished...
for ( $row = 0; $row < $rows; ++$row )
if ( ! $isDone[ $row ] )
{
$firstTerm = $matrix[ $row ][ $column ];
// Back substitution.
for ( $subColumn = $column; $subColumn <= $rows; ++$subColumn )
{
$accumulator =
bcmul( $firstTerm, $matrix[ $activeRow ][ $subColumn ] );
$matrix[ $row ][ $subColumn ] =
bcsub( $matrix[ $row ][ $subColumn ], $accumulator );
}
}
// Save this coefficient for the return.
$coefficients[ $column ] = $matrix[ $activeRow ][ $rows ];
}
// Coefficients are stored backward, so sort them.
ksort( $coefficients );
// Return the coefficients.
return $coefficients;
} // getCoefficients
static public function interpolate( $coefficients, $x )
{
$numberOfCoefficient = count( $coefficients );
$y = 0;
for ( $coefficentIndex = 0; $coefficentIndex < $numberOfCoefficient; ++$coefficentIndex )
{
// y += coefficients[ coefficentIndex ] * x^coefficentIndex
$y =
bcadd
(
$y,
bcmul
(
$coefficients[ $coefficentIndex ],
bcpow( $x, $coefficentIndex )
)
);
}
return floatval( $y );
} // interpolate
} // Class
テスト実行
データは、配列で指定する。 中の array(1,0),・・・が、array(Xの値,Yの値（目的変数）)の順である。
require_once( 'inc_PolynomialRegression.php' ); //class読込
//テスト配列定義
$data = array (
array(1,0),
array(2,4),
array(3,3),
array(4,2),
array(5,5),
array(6,3),
array(7,8)
);
//print_r ( $data ) ;
// Precision digits in BC math.
bcscale( 10 );
// Start a regression class of order 2--linear regression.
$PolynomialRegression = new PolynomialRegression( 2 ); //変数の数、次数+1
// Add all the data to the regression analysis.
foreach ( $data as $dataPoint )
$PolynomialRegression->addData( $dataPoint[ 0 ], $dataPoint[ 1 ] );
// Get coefficients for the polynomial.
$coefficients = $PolynomialRegression->getCoefficients();
//
// Get average of Y-data.
//
$Y_Average = 0.0;
foreach ( $data as $dataPoint )
$Y_Average += $dataPoint[ 1 ];
$Y_Average /= count( $data );
//
// Calculate R Squared.
//
$Y_MeanSum = 0.0;
$Y_ErrorSum = 0.0;
foreach ( $data as $dataPoint )
{
$x = $dataPoint[ 0 ];
$y = $dataPoint[ 1 ];
$error = $y;
$error -= $PolynomialRegression->interpolate( $coefficients, $x );
$Y_ErrorSum += $error * $error;
$error = $y;
$error -= $Y_Average;
$Y_MeanSum += $error * $error;
}
$R_Squared = 1.0 - ( $Y_ErrorSum / $Y_MeanSum );
// Print slope and intercept of linear regression.
// 四捨五入 4桁目
$para_a = round( $coefficients[ 1 ], 4 );
$para_b = round( $coefficients[ 0 ], 4 );
$para_R = round( $R_Squared ,4 );
//結果出力
print "$para_a,$para_b,$para_R" ;
出力結果
0.8571,0.1429,0.5455
よって
ｙ＝0.8571 ｘ + 0.1429
重決定 R2 は 0.5455
※X,Yの入れ替えに注意
エクセルとの比較：OK
二次関数で近似するときは
$PolynomialRegression = new PolynomialRegression( 3 );//2→3にする $para_a = round( $coefficients[ 2 ], 4 );//増やす $para_b = round( $coefficients[ 1 ], 4 ); $para_c = round( $coefficients[ 0 ], 4 ); $para_R = round( $R_Squared ,4 );
出力：0.0952,0.0952,1.2857,0.5657
エクセルと比較：OK
/*=========================================================================*/ /* Name: PolynomialRegression.php */ /* Uses: Calculates and returns coefficients for polynomial regression. */ /* Date: 06/01/2009 */ /* Author: Andrew Que (http://www.DrQue.net/) */ /* Revisions: */ /* 0.8 - 06/01/2009- QUE - Creation. */ /* 0.9 - 06/14/2012- QUE - */ /* + Bug fix: removed notice causes by uninitialized variable. */ /* + Converted naming convention. */ /* + Fix spelling errors (or the ones I found). */ /* + Changed to row-echelon method for solving matrix which is much */ /* faster than the determinant method. */ /* 0.91 - 05/17/2013- QUE - */ /* = Changed name to Polynonial regression as this is more fitting to */ /* to the function. */ /* 0.92 - 12/28/2013- QUE - */ /* + Added forced offset. */ /* 1.00 - 12/29/2013 - QUE - */ /* + Forced offset changed to allow any term to be forced. */ /* Unit complete. Correlation coefficient (r-squared) has been */ /* implemented externally in the demos. */ /* 1.1 - 2014/05/05 - QUE - */ /* + 'interpolate' is now static as it does not need an instance to */ /* operate. Useful if coefficients have been calculated elsewhere. */ /* - Deprecated 'setDegree' function. This is the wrong terminology for */ /* what the function does. It actually sets the number of */ /* coefficients for the polynomial. The degree of the polynomial is */ /* the number of coefficients less one. Made the identical function */ /* 'setNumberOfCoefficient' to replace it. */ /* + Added getter functions for anything that has a set function. */ /* */ /* This project is maintained at: */ /* http://PolynomialRegression.drque.net/ */ /* */ /* ----------------------------------------------------------------------- */ /* */ /* Polynomial regression class. */ /* Copyright (C) 2009, 2012-2014 Andrew Que */ /* */ /* This program is free software: you can redistribute it and/or modify */ /* it under the terms of the GNU General Public License as published by */ /* the Free Software Foundation, either version 3 of the License, or */ /* (at your option) any later version. */ /* */ /* This program is distributed in the hope that it will be useful, */ /* but WITHOUT ANY WARRANTY; without even the implied warranty of */ /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */ /* GNU General Public License for more details. */ /* */ /* You should have received a copy of the GNU General Public License */ /* along with this program. If not, see
* Used for calculating polynomial regression coefficients. Useful for
* linear and non-linear regression, and polynomial curve fitting.
*
* @package PolynomialRegression
* @author Andrew Que ({@link http://www.DrQue.net/})
* @link http://PolynomialRegression.drque.net/ Project home page.
* @copyright Copyright (c) 2009, 2012-2014, Andrew Que
* @license http://opensource.org/licenses/gpl-license.php GNU Public License
* @version 1.1
*/
/**
* Used for calculating polynomial regression coefficients and interpolation using
* those coefficients. Useful for linear and non-linear regression, and polynomial
* curve fitting.
*
* Note: Requires BC math to be compiled into PHP. Higher-degree polynomials end up
* with very large/small numbers, requiring an arbitrary precision arithmetic. Make sure
* to set "bcscale" as coefficients will likely have decimal values.
*
* Quick example of using this unit to calculate linear regression (1st degree polynomial):
*
*
* $regression = new PolynomialRegression( 2 );
* // ...
* $regression->addData( $x, $y );
* // ...
* $coefficients = $regression->getCoefficients();
* // ...
* $y = $regression->interpolate( $coefficients, $x );
*
*
*
* @package PolynomialRegression
* @link http://PolynomialRegression.drque.net/ Project home page.
0 件のコメント:
コメントを投稿