Location: PHPKode > scripts > Numerical > numerical/class.Numerical.php
<?php

//
// Edit History:
//
//  $Author: munroe $
//  $Date: 2006/01/27 15:48:36 $
//
//  Dick Munroe (hide@address.com) 04-Jan-2006
//      Initial version created.
//
//  Dick Munroe (hide@address.com) 04-Jan-2006
//      Add random variable generation using the Box-Mueller transformations.
//
//  Dick Munroe (hide@address.com) 05-Jan-2006
//      Add a mean, median, and mode functions.
//
//  Dick Munroe (hide@address.com) 06-Jan-2006
//      Add standard deviation and variance functions.
//
//  Dick Munroe (hide@address.com) 19-Jan-2006
//      Improve the numerical accuracy and performance of the variance and
//      standard deviation functions (suggested by Martin Weis).
//
//  Dick Munroe (hide@address.com) 22-Jan-2006
//      Swap the integration limits internally to the method if they are
//      out of order.
//
//  Dick Munroe (hide@address.com) 25-Jan-2006
//      Add an equation solver using steepest descent and another using
//      bisection.
//
//  Dick Munroe (hide@address.com) 26-Jan-2006
//      Add a limiting feature to the bisection equation solver.
//
//  Dick Munroe (hide@address.com) 27-Jan-2006
//      Add conversion to rational numbers.
//
//  Dick Munroe (hide@address.com) 28-Jan-2006
//      Add sieve of Erastaphenes.
//      Add prime factoring.
//      Use prime factoring to reduce rational numbers to lowest form.
//

/**
 * @author Dick Munroe <hide@address.com>
 * @copyright copyright @ 2006 by Dick Munroe, Cottage Software Works, Inc.
 * @license http://www.csworks.com/publications/ModifiedNetBSD.html
 * @version 1.0.8
 * @package Numerical
 * @example ./example.php
 */

/*
 * Provide a number of numerical analysis and statistical functions.
 */

class Numerical
{
    /**
     * Factor a number into it's primes.
     *
     * @access public
     * @param integer $n the number to be factored.
     * @return reference to array of integer prime factors are the keys, the values are the number
     *                          times the factor occurs.
     */

    function &factor($n)
    {
        $x = ceil(sqrt(abs($n))) ;

        $thePrimes = Numerical::sieve($x) ;

        if ($thePrimes == NULL)
        {
            return array(abs($n) => 1) ;
        }

        foreach ($thePrimes as $aPrime)
        {
            do
            {
                if (($n % $aPrime) == 0)
                {
                    $theFactors[$aPrime]++ ;
                    $n = $n / $aPrime ;
                    if ($n == 1)
                    {
                        return $theFactors ;
                    }
                }
                else
                {
                    break ;
                }
            } while (true) ;
        }

        $theFactors[$n]++ ;

        return $theFactors ;
    }

    /**
     * Generate a Gaussian normal distribution with a specified mean and standard
     * deviation.  The default distribution generated is the standard
     * normal distribution with a mean of 0.0 and a standard deviation of
     * 1.0.
     *
     * @desc Draw a value from a normal distribution.
     * @access public
     * @param float $x The parameter of the function.
     * @param float $theMean The mean of the distribution, by default 0.0.
     * @param float $theStandardDeviation The standard deviation of the distribution, by default, 1.0.
     * @return float the probability of the x occuring in the distribution.
     */

    function gaussian($x, $theMean = 0.0, $theStandardDeviation = 1.0)
    {
        $part1 = (1/($theStandardDeviation * sqrt(2 * M_PI))) ;
        $part2 = exp(- pow(($x - $theMean), 2) / (2 * pow($theStandardDeviation, 2))) ;
        return $part1 * $part2 ;
    }

    /*
     * Caluclate the integral of a function over a specified range.
     *
     * Uses the midpoint method of calculation for integral.
     * The function passed in must take one parameter and return the
     * value of the function for that parameter.
     *
     * @access public
     * @param function $theFunction the function to be integrated.
     * @param float $theLowLimit The lower limit of the integral
     * @param float $theHighLimit The high limit of the integral
     * @param integer $theNumberOfSteps The number of steps take to evaluate the integral, defaults to 100.
     * @return float The value of the integral for the specified range.
     */

    function integrate($theFunction, $theLowLimit, $theHighLimit, $theNumberOfSteps = 100)
    {
        if ($theLowLimit > $theHighLimit)
        {
            /*
             * The limits are out of order, swap them.
             */

            $xxx = $theLowLimit ;
            $theLowLimit = $theHighLimit ;
            $theHighLimit = $xxx ;
        }

        $theDelta = ($theHighLimit - $theLowLimit) / $theNumberOfSteps ;
        $theMidpointDelta = $theDelta / 2 ;

        $theArea = 0.0 ;

        for ($i = 0; $i < $theNumberOfSteps; $i++)
        {
            $theValue = $theLowLimit + (($i * $theDelta) + $theMidpointDelta) ;
            $theValue = $theFunction($theValue) ;
            $theArea += $theValue * $theDelta ;
        }

        return $theArea ;
    }

    /*
     * @desc Calculate the mean of an array of samples.
     * @access public
     * @param array $theSamples An array of samples.
     * @return float the mean of the samples.
     */

    function mean($theSamples)
    {
        return array_sum($theSamples) / count($theSamples) ;
    }

    /*
     * @desc Determine the median of a sample.
     * @access public
     * @param array $theSamples An array of samples.
     * @return float the mean of the samples.
     */

    function median($theSamples)
    {
        sort($theSamples) ;

        if ((count($theSamples) % 2) == 0)
        {
            $i = count($theSamples) / 2 ;
            return ($theSamples[$i - 1] + $theSamples[$i]) / 2 ;
        }
        else
        {
            return $theSamples[(int) (count($theSamples) / 2)] ;
        }
    }

    /*
     * @desc Determine the mode[s] of a sample.
     * @access public
     * @param array $theSamples An array of samples.
     * @return mixed the mode[s] of the sample.
     */

    function mode($theSamples)
    {
        $theCounts = array() ;

        $theCount = count($theSamples) ;

        for ($i = 0; $i < $theCount; $i++)
        {
            $theCounts[(string)$theSamples[$i]]++ ;
        }

        arsort($theCounts) ;

        $theModes = array() ;

        foreach ($theCounts as $key => $count)
        {
            if (count($theModes) == 0)
            {
                $theModeCount = $count ;
            }

            if ($theModeCount == $count)
            {
                $theModes[] = $key ;
            }
            else
            {
                break ;
            }
        }

        return $theModes ;
    }

    /*
     * Note that this algorithm can be unstable for uniformly distributed random numbers
     * close to 0.0.
     *
     * @desc Generate gaussian distributed random numbers using the Box-Mueller basic transform.
     * @param float $theMean The mean of the distribution.
     * @param float $theStandardDeviation The standard deviation of the distribution.
     * @return float A random number.
     */

    function randomGaussianBoxMuellerBasic($theMean = 0.0, $theStandardDeviation = 1.0)
    {
        static $useLast = FALSE ;
        static $y2 ;

        if ($useLast)
        {
            $useLast = FALSE ;
            $y1 = $y2 ;
        }
        else
        {
            /*
             * Generate a uniformly distributed random number in the range
             * (0, 1].
             */

            $theRange = 10000000 ;

            $x1 = mt_rand(1, $theRange) / $theRange ;
            $x2 = mt_rand(1, $theRange) / $theRange ;

            /*
             * Convert this to a pair of gaussian distributed random numbers.
             * The 2nd is returned on the next call.
             */

            $x = sqrt( -2 * log($x1) ) ;

            $y1 = $x * sin( 2 * M_PI * $x2 ) ;
            $y2 = $x * cos( 2 * M_PI * $x2 ) ;

            $useLast = TRUE ;
        }

        return $theMean + $y1 * $theStandardDeviation ;
    }

    /*
     * @desc Generate gaussian distributed random numbers using the Box-Mueller Polar transform.
     * @param float $theMean The mean of the distribution.
     * @param float $theStandardDeviation The standard deviation of the distribution.
     * @return float A random number.
     */

    function randomGaussianBoxMuellerPolar($theMean = 0.0, $theStandardDeviation = 1.0)
    {
        static $useLast = FALSE ;
        static $y2 ;

        if ($useLast)
        {
            $useLast = FALSE ;
            $y1 = $y2 ;
        }
        else
        {
            do
            {
                /*
                 * Generate a pair of uniformly distributed random numbers
                 * in the range [-1..1].
                 */

                $theRange = 10000000 ;

                $x1 = mt_rand(- $theRange, $theRange) / $theRange ;
                $x2 = mt_rand(- $theRange, $theRange) / $theRange ;

                $w = $x1 * $x1 + $x2 * $x2 ;
            } while ($w >= 1.0) ;

            $w = sqrt ( (-2 * log( $w ) ) / $w ) ;
            $y1 = $x1 * $w ;
            $y2 = $x2 * $w ;

            $useLast = TRUE ;
        }

        return $theMean + $y1 * $theStandardDeviation ;
    }

    /*
     * @desc Determine the standard deviation of a sample.
     * @access public
     * @param array $theSamples An array of samples.
     * @return float the standard deviation.
     */

    function standardDeviation($theSamples)
    {
        return sqrt(Numerical::variance($theSamples)) ;
    }

    /*
     * Calculate the rational number version of a floating point number.
     *
     * Note that the rational number returned is guaranteed to be in
     * simplest form, i.e., if the number discovered is of the form 4/6,
     * it will be reduced to 2/3 prior to return..
     *
     * A very cool algorithm (and, as near as I can tell, unique).
     *
     *  1. Take a guess at a denominator.
     *  2. Take a guess at a numerator.
     *  3. Check the numerator / denominator to see how close we came to the
     *     input vaue.
     *  4. Calculate the "remaining" portion of the input number not
     *     covered by the numerator / denominator.
     *  5. Take a guess at a denominator for the remainder.
     *  6. Take a guess at the numerator for the input number, given the
     *     new gess for the denominator.
     *  7. Go to 3.
     *
     * @param float $x The floating point number to be converted.
     * @param float $theEpsilon How close the fraction should be to x.
     * @param float $theLimit How big the numberator or denominator can get.
     * @return array [0] = numerator [1] = denominator.
     */

    function rational($x, $theEpsilon = 1.0e-06)
    {
        /*
         * Strip off the sign, the algorithm only works on positive
         * numbers.  We'll put the sign back on the way out.
         */

        if ($x < 0.0)
        {
            $theSign = -1 ;
            $x = -$x ;
        }
        else
        {
            $theSign = 1 ;
        }

        /*
         * Strip off the integer portion, we'll add that back in as we return.
         */

        $theInteger = floor($x) ;
        $x = $x - $theInteger ;

        /*
         * Catch the case of a real number with no fractional part.
         */

        if ($x < $theEpsilon)
        {
            return array((int)($theSign * $theInteger), 1) ;
        }

        /*
         * Take an inital guess at the fractional part of the rational number.
         */

        $d = round(1/$x, 0) ;
        $n = round($x * $d, 0) ;

        do
        {
            /*
             * If our current ratio is "close enough", we're done.
             */

            $delta = abs($x - ($n / $d)) ;

            if ($delta < $theEpsilon)
            {
                $theFactors = Numerical::factor($n) ;

                foreach ($theFactors as $aFactor => $theRepeat)
                {
                    for ($i = 0; $i < $theRepeat; $i++)
                    {
                        if (($d % $aFactor) == 0)
                        {
                            $n = $n / $aFactor ;
                            $d = $d / $aFactor ;
                        }
                        else
                        {
                            breqk ;
                        }
                    }
                }

                $n = $theInteger * $d + $n ;

                /*
                 * return error if the numerator overflows.
                 */

                if ($n < 0)
                {
                    return sqrt(-1) ;
                }

                return array((int)($theSign * $n), (int)$d) ;
            }

            /*
             * Figure out the denominator of the fraction representing
             * the delta.
             */

            $d = round(1/$delta, 0) ;

            /*
             * Check to see if we had an integer overflow in calculating
             * the denominator.  If we did, then bail with Nan.
             */

            if ($d < 0)
            {
                return sqrt(-1) ;
            }

            /*
             * The new guess for the numerator is however many pieces of
             * the denominator are necessary.
             */

            $n = round($x * $d, 0) ;
        } while (true) ;
    }

    /**
     * Sieve of Erostophenes.
     *
     * @access public
     * @param integer $n the maximum number for which the sieve is to be generated.
     * @return array of integer All primes in the range 2..n.
     */

    function sieve($n)
    {
        if ($n < 2)
        {
            return NULL ;
        }

        /*
         * Populate the candidate values for primes.
         */

        $theSieve[2] = 2 ;

        for ($i = 3; $i <= $n; $i = $i + 2)
        {
            $theSieve[$i] = $i ;
        }

        for ($i = 3; $i <= $n; $i = $i + 2)
        {
            if (isset($theSieve[$i]))
            {
                for ($j = $i + $i; $j <= $n; $j = $j + $i)
                {
                    unset($theSeive[$j]) ;
                }
            }
        }

        return array_keys($theSieve) ;
    }

    /*
     * Solve for roots using the Bisection algorithm.
     *
     * The value of the function must change sign from the left to the
     * right hand side.  This method is slow, but certain.
     *
     * @access public
     * @param function $theFunction the function to be solved.
     * @param float $theLHS The left hand side of the bracket.
     * @param float $theRHS The right hand side of the bracket.
     * @param float $theEpsilon How close does the solution have to be to zero to be a solution.
     * @param float $theSSEpsilon How small does the solution space have to be before failure is declared.
     * @return float The value of the input that solves the function.
     */

    function solveBisection($theFunction, $theLHS, $theRHS, $theEpsilon = 1.0e-06, $theSSEpsilon = 1.0e-9)
    {
        while (true)
        {
            /*
             * If the solution space has gotten very small and no
             * solution is within an epsilon of 0.0, then we declare
             * failure.
             */

            if (($theRHS - $theLHS) < $theSSEpsilon)
            {
                return sqrt(-1) ;
            }

            $x = ($theLHS + $theRHS) / 2.0 ;

            /*
             * Calculate the value of the function.
             */

            $y = $theFunction($x) ;

            /*
             * Since this always looks for 0, if we're within epsilon of any 0,
             * quit.
             */

            if (abs($y) < $theEpsilon)
            {
                return $x ;
            }

            /*
             * Switch the appropriate side of the bracket to use
             * the guessed $x value for a root.
             */

            if ($y < 0)
            {
                $theLHS = $x ;
            }
            else
            {
                $theRHS = $x ;
            }
        }
    }

    /*
     * Solve for roots using the steepest descent algorithm.
     *
     * Can iterate a lot for "flat" curves.
     *
     * @access public
     * @param function $theFunction the function to be solved.
     * @param float $theGuess The guess at a solution.
     * @param float $theEpsilon Accuracy of the solution.
     * @param float $theDerivativeDelta Delta over which the derivative is calculated.
     * @return float The value of the input that solves the function.
     */

    function solveSteepestDescent($theFunction, $theGuess, $theEpsilon = 1.0e-06, $theDerivativeDelta = .001,
                                  $theIterationLimit = 20)
    {
        while (true)
        {
            /*
             * Calculate the value of the function.
             */

            $y = $theFunction($theGuess) ;

            /*
             * Since this always looks for 0, if we're within epsilon of any 0,
             * quit.
             */

            if (abs($y) < $theEpsilon)
            {
                return $theGuess ;
            }

            /*
             * There are some function which can take a LONG
             * time to converge, so we limit the number of
             * iterations and return NaN if we get into this
             * situation.
             */

            if ($theIteration++ > $theIterationLimit)
            {
                return @sqrt(-1) ;
            }

            /*
             * The derivative is the slope "at" the point in question ($theGuess).
             * We define the derivative delta as that distance which is considered
             * for all practical purposes, 0.
             */

            $theDerivative = ($theFunction($theGuess + $theDerivativeDelta) - $y) / $theDerivativeDelta ;

            /*
             * The new guess is calculate by assuming that the
             * equation being solved is a line and moving to where
             * the solution would be.  The algegra is pretty easy
             * keeping in mind that the derivative is the slope of
             * the curve at the point $y.  In effect this subtracts
             * the solution for the linear approximation of the
             * curve at the point $theGuess from $theGuess, moving
             * in the direction of the solution.
             */

            $theGuess = $theGuess - ($y / $theDerivative) ;
        }
    }

    /*
     * @desc Determine the variance of a sample.
     * @access public
     * @param array $theSamples An array of samples.
     * @return float the variance.
     */

    function variance($theSamples)
    {
        $theMean = Numerical::mean($theSamples) ;
        $theSum = 0.0 ;
        foreach ($theSamples as $theValue)
        {
            $theSum += pow(($theValue - $theMean), 2) ;
        }
        return ($theSum/(count($theSamples) - 1)) ;
     }
}

?>
Return current item: Numerical