#b-navbar { height:0px; display:none; visibility:hidden; }

ページ

2014年6月22日日曜日

PHP:単回帰分析(係数と重決定R2値)(1次関数~n次関数)

まず1次関数での近似を行った。

出典: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

よって
y=0.8571 x + 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 . */ /* */ /* ----------------------------------------------------------------------- */ /* */ /* (C) Copyright 2009, 2012-2014 */ /* Andrew Que */ /* > */ /*=========================================================================*/ /** * Polynomial regression. * *

* 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 件のコメント:

コメントを投稿