Gauss-Newton class
A PHP non-linear regression class
Introduction
The Gauss-Newton method is a non-linear curve fitting algorithm. It uses the minimization of least-squares in order to determine coefficients that can be used to in a function that describes data. This method can be used on any differentiable function on a data set of arbitrary size.
The process involves making an initial guess at the function coefficients, and then using the Gauss-Newton method to refine them. As an iterative process based on an initial guess means the algorithm may never converge. It may also oscillate wildly or diverge if the initial guess is too far from the actual coefficients. This often is not an issue as the information about the function being modeled and the data coming in are often fairly define. This helps assemble an initial guess that should prevent divergence.
This PHP class is setup in two parts. The first part is the Gauss-Newton function itself, which is an abstract class. It defines one function that can be used to refine the coefficients. The additional classes are implementations of various mathematical functions that can be used for fitting.
Additional functions are fairly easy to create. Methods are define that can return a value for a partial derivative with respect to one of the coefficients, as well as a method for the function itself.
There are two versions of this library. One uses conventional floating-point values, and the other uses BC Math. They are identical in methods with the only difference being the number system used for calculations. Naturally, the BC Math implementation is significantly slower but can deal with larger numbers and large sample sizes.
Unmaintained
This library hasn't proved to be very useful and isn't maintained. You can still download and play with it, but the library isn't recommended.
Functions currently implemented:
- f( x ) = a eb x
- f( x ) = a eb x+ c
- f( x ) = a eb x + c + d (Not useful—never converges)
- f( x ) = a ( 1 - eb x )
- f( x ) = a ( 1 - eb x ) + c
- f( x ) = a ( 1 - eb x + c ) + d (Not useful—never converges)
- f( x ) = a eb x sin( c x + d )
- f( x ) = a eb x ( cos( c x ) + sin( c x ) )
- f( x ) = a eb x ( cos( c x + d ) + sin( c x + d ) )
- f( x ) = a xb
- f( x ) = a xb + c
- f( x ) = a nb x
- f( x ) = a nb x + c
- f( x ) = a bc x (Not useful—never converges)
- f( x ) = a bc x + d (Not useful—never converges)
- f( x ) = a bc x + d + g (Not useful—never converges)
Manual
Documentation is available online, generated by from the source code using phpDocumentor.
Download
Current versions
Version 1.1, Released June 2, 2013.
This version has some speed improvements to the BC Math trigonometric function.
BZip 2 (Linux)
Download release 1.1MD5: 480ad868484673c4efa357373eee1c1f
SHA1: 2421d948667444198081fdb0cee5fde1df8f65bc
Zip (Others)
Download release 1.1MD5: d5d07053ec0c98f4464a5a3043a271a3
SHA1: 8c5e08333861b68f6d09c4e5648fa02fe8485d05
Previous versions
Version 1.0, Released May 27, 2013.
BZip 2 (Linux)
Download release 1.0MD5: 1bb1f06810cf1a651e31042feb141e2e
SHA1: 1091ba8b8546b47653f3158ca8064ea28e2382fa
Zip (Others)
Download release 1.0MD5: 1bdfa8cceb92d3beb2ea56389a1215f1
SHA1: 62b13a77eb7f8488fcbcb17d639b1dacd3b894e1
User comments
From Andrew Que (http://www.DrQue.net/)
Middleton, WI, USA
May 27, 2013 at 5:50 PM
Examples
Basic Exponential Regression.
Here a set of data points that perfectly matches an exponential curve. The coefficients for this data are unknown, and the Gauss-Newton method is used to calculate them.
<?php
// Load the least-square regression class.
require_once( 'RootDirectory.inc.php' );
require_once( $RootDirectory . 'Includes/GaussNewton/ExponentialRegression.php' );
require_once( $RootDirectory . 'Includes/plot.php' );
// Maximum number of refinement iterations.
$MAX_ITERATIONS = 20;
// Data to curve fit.
// True function is f( x ) = 2.5 * exp( -1.2 * x ).
$data =
array
( // x y
array( 0, 2.50000000000000 ),
array( 1, 0.75298552978051 ),
array( 2, 0.22679488322353 ),
array( 3, 0.06830930611823 ),
array( 4, 0.02057436762255 ),
array( 5, 0.00619688044167 ),
);
// Initial guess at coefficients.
$coefficients = array( 1, -1 );
$lastCoefficients = array( 0, 0 );
// Create an instance of the regression to use.
$regression = new ExponentialRegression();
// Counter for the number of iterations.
$iteration = 0;
// Refine coefficients.
while ( ( $coefficients != $lastCoefficients )
&& ( $iteration < $MAX_ITERATIONS ) )
{
$lastCoefficients = $coefficients;
// Print the results.
echo "Iteration $iteration: f( x ) = $coefficients[0] * exp( $coefficients[1] * x )<br />";
// Run Gauss-Newton method with current coefficients.
$coefficients =
$regression->refineCoefficients
(
array_column( $data, 0 ),
array_column( $data, 1 ),
$coefficients
);
++$iteration;
}
?>
The output of this example:
Iteration 1: f( x ) = 2.4981198672862 * exp( -1.4286924972401 * x )
Iteration 2: f( x ) = 2.499368425943 * exp( -1.1633050892343 * x )
Iteration 3: f( x ) = 2.4999624788228 * exp( -1.1990654527662 * x )
Iteration 4: f( x ) = 2.4999999782808 * exp( -1.1999994034099 * x )
Iteration 5: f( x ) = 2.5 * exp( -1.1999999999998 * x )
Iteration 6: f( x ) = 2.5 * exp( -1.2 * x )
Note that the code runs the refine function until the coefficients no longer change. This happens when the function converges to the correct answer, or when the floating-point numbers run out of precision. Because the initial guesses are fairly close to the true coefficients, it only take 6 iterations to obtain the exact values.
There are cases when the result will not converge so quickly, or may never converge. To cover such situations, a maximum number of iterations should be enforced.
Regression Graphed.
This time the data has noise induced error. The function plot that is used to draw the graph can be viewed here. It uses XY Plot to draw the chart.
<?php
require_once( 'RootDirectory.inc.php' );
require_once( $RootDirectory . 'Includes/GaussNewton/ExponentialRegression.php' );
require_once( $RootDirectory . 'Includes/plot.php' );
// Maximum number of refinement iterations.
$MAX_ITERATIONS = 20;
// Data to curve fit.
// True function is f( x ) = 5.2 * exp( -0.13 * x ) with 0.2 standard
// deviation noise.
$data =
array
( // x y
array( 0, 5.5144807349 ), array( 11, 1.1326357272 ),
array( 1, 4.7759860442 ), array( 12, 1.3408739643 ),
array( 2, 3.8420020059 ), array( 13, 1.057980605 ),
array( 3, 3.5628540333 ), array( 14, 1.0143909663 ),
array( 4, 3.0456179633 ), array( 15, 0.8429361305 ),
array( 5, 3.0429401144 ), array( 16, 0.3357937422 ),
array( 6, 2.1533753396 ), array( 17, 0.4535925247 ),
array( 7, 2.4085397215 ), array( 18, 0.6251204346 ),
array( 8, 1.6196966715 ), array( 19, 0.4746199103 ),
array( 9, 1.5311146956 ), array( 20, 0.3118236578 ),
array( 10, 1.3331911728 ),
);
// Initial guess at coefficients.
$coefficients = array( 1, -0.05 );
$lastCoefficients = array( 0, 0 );
// Create an instance of the regression to use.
$regression = new ExponentialRegression();
// Counter for the number of iterations.
$iteration = 0;
// Refine coefficients.
while ( ( $coefficients != $lastCoefficients )
&& ( $iteration < $MAX_ITERATIONS ) )
{
$lastCoefficients = $coefficients;
// Run Gauss-Newton method with current coefficients.
$coefficients =
$regression->refineCoefficients
(
array_column( $data, 0 ),
array_column( $data, 1 ),
$coefficients
);
++$iteration;
}
// Create a string with the derived function.
$functionString =
"f( x ) = "
. round( $coefficients[ 0 ], 5 )
. " * exp( "
. round( $coefficients[ 1 ], 5 )
. " * x ), Iterations: $iteration";
// Plot this data.
plot
(
$data,
$regression,
$coefficients,
$functionString
);
?>
The output of this example:
The calculated coefficients do not exactly match the coefficients used to generate the data because of the noise that has been added. Nonetheless the curve fits the data well.
Uneven Interval Regression.
The spacing between x values does not need to be regular intervals. This code shows a concentration of samples near the beginning of the function, but few at the end. The algorithm is still able to deduce the coefficients.
This can be useful for situations where data may not be accurate unless the amplitude is large.
<?php
require_once( 'RootDirectory.inc.php' );
require_once( $RootDirectory . 'Includes/GaussNewton/SineRegression.php' );
require_once( $RootDirectory . 'Includes/plot.php' );
// Maximum number of refinement iterations.
$MAX_ITERATIONS = 20;
// Data to curve fit.
// True function is
// f( x ) = 2.2 * exp( -0.045 * x ) * ( cos( 0.16 * pi * x ) + sin( 0.16 * pi * x ) )
$data =
array
( // x y
array( 100.0000000000000, 0.0244397923841 ), array( 14.6214911951032, 1.5481465627215 ),
array( 65.0514997831991, 0.1464156478272 ), array( 12.9818655252878, 1.4850366650802 ),
array( 52.2878745280338, 0.2763956423459 ), array( 11.4573994178928, 0.4800179585686 ),
array( 44.3028323846582, -0.3703441442617 ), array( 10.0329725273209, -0.8699313691027 ),
array( 38.4775539310863, 0.5271225713246 ), array( 8.6962598649587, -1.8995252010157 ),
array( 33.8890352633040, -0.5803394944407 ), array( 7.4370825640463, -2.1868249308332 ),
array( 30.1029995663981, -0.1664584025433 ), array( 6.2469368304150, -1.6583168013534 ),
array( 26.8801001050522, 0.9163370418915 ), array( 5.1186454354779, -0.5313631966176 ),
array( 25.4319153082864, 0.8349306763142 ), array( 4.5757490560675, 0.1421691299848 ),
array( 22.7965977824862, -0.3531355887698 ), array( 3.5290537142854, 1.4597816617840 ),
array( 20.4467696486751, -1.2369479914185 ), array( 2.5304996677544, 2.4541849094176 ),
array( 18.3265772210207, -0.7389345507815 ), array( 1.5758525723033, 2.8982150892254 ),
array( 16.3951071032141, 0.5770189809001 ), array( 0.6614132866878, 2.7155109372628 ),
);
// Initial guess at coefficients.
$coefficients = array( 1, -0.01, .4 );
$lastCoefficients = array( 0, 0, 0 );
// Create an instance of the regression to use.
$regression = new DampedSineRegression2();
// Counter for the number of iterations.
$iteration = 0;
// Refine coefficients.
while ( ( $coefficients != $lastCoefficients )
&& ( $iteration < $MAX_ITERATIONS ) )
{
$lastCoefficients = $coefficients;
// Run Gauss-Newton method with current coefficients.
$coefficients =
$regression->refineCoefficients
(
array_column( $data, 0 ),
array_column( $data, 1 ),
$coefficients
);
++$iteration;
}
// Create a string with the derived function.
$functionString =
"f( x ) = "
. round( $coefficients[ 0 ], 3 )
. " * exp( "
. round( $coefficients[ 1 ], 3 )
. " * x ) * ( sin( "
. round( $coefficients[ 2 ], 3 )
. " * x ) + cos( "
. round( $coefficients[ 2 ], 3 )
. " * x ) )";
// Plot this data.
plot
(
$data,
$regression,
$coefficients,
$functionString
);
?>
The output of this example:
Custom Curve.
To implement a new method of regression requires creating a child class of GaussNewton and defining the abstract functions. The following example demonstrates an implementation of inverse square with a scale and offset coefficient.
<?php
require_once( 'RootDirectory.inc.php' );
require_once( $RootDirectory . 'Includes/GaussNewton/GaussNewtonRegression.php' );
require_once( $RootDirectory . 'Includes/plot.php' );
class InverseSquareRegression extends GaussNewtonRegression
{
//
// Partial differential.
//
protected function partialDifferential( $x, $coefficientIndex, $coefficients )
{
$result = 0;
switch ( $coefficientIndex )
{
// Partial derivative with respect to first coefficient.
case 0:
{
$result = 1 / pow( $x, 2 );
break;
}
// Partial derivative with respect to second coefficient.
case 1:
{
$result = 1;
break;
}
// If this function is called with some other index, force an assertion.
default:
{
assert( false );
break;
}
}
return $result;
}
//
// Number of coefficients.
//
protected function getNumberOfCoefficients()
{
return 2;
}
//
// Evaluate function at x.
//
public function getFunction( $x, $coefficients )
{
return $coefficients[ 0 ] / pow( $x, 2 ) + $coefficients[ 1 ];
}
}
// Maximum number of refinement iterations.
$MAX_ITERATIONS = 20;
// Data to curve fit.
// True function is
// f( x ) = 2.6 / x^2 + 1.2
$data =
array
( // x y
array( 1, 3.8 ),
array( 2, 1.85 ),
array( 3, 1.4888888889 ),
array( 4, 1.3625 ),
array( 5, 1.304 ),
);
// Initial guess at coefficients.
$coefficients = array( 1, 0 );
$lastCoefficients = array( 0, 0 );
// Create an instance of the regression to use.
$regression = new InverseSquareRegression();
// Counter for the number of iterations.
$iteration = 0;
// Refine coefficients.
while ( ( $coefficients != $lastCoefficients )
&& ( $iteration < $MAX_ITERATIONS ) )
{
$lastCoefficients = $coefficients;
// Run Gauss-Newton method with current coefficients.
$coefficients =
$regression->refineCoefficients
(
array_column( $data, 0 ),
array_column( $data, 1 ),
$coefficients
);
++$iteration;
}
// Create a string with the derived function.
echo "f( x ) = "
. round( $coefficients[ 0 ], 3 )
. " / x^2 + "
. round( $coefficients[ 1 ], 3 );
?>
The output of this example:
The function a / x2 + b is define with getFunction. The coefficient a is coefficient[ 0 ], and b is coefficient[ 1 ]. There are two coefficients, and therefore two partial derivatives. The partial derivative with respect to coefficient a is 1 / x2, and with respect to coefficient b simply 1 for all x. These are defined in partialDifferential.
Note that the inverse square has already been implemented with FixedPowerRegression2. For this example, n is would be -2. The second coefficient to scale x would always be 1.
License
This software is free, open-source software released under the GNU license.
Author
The Gauss-Newton class is written and maintained by Andrew Que. To get in touch with Andrew Que, visit his contact page