8889841cStudentT.php 0000644 00000007720 15053602601 0007033 0 ustar 00 getMessage();
}
if (($value < 0) || ($degrees < 1) || ($tails < 1) || ($tails > 2)) {
return Functions::NAN();
}
return self::calculateDistribution($value, $degrees, $tails);
}
/**
* TINV.
*
* Returns the one-tailed probability of the chi-squared distribution.
*
* @param mixed $probability Float probability for the function
* @param mixed $degrees Integer value for degrees of freedom
*
* @return float|string The result, or a string containing an error
*/
public static function inverse($probability, $degrees)
{
$probability = Functions::flattenSingleValue($probability);
$degrees = Functions::flattenSingleValue($degrees);
try {
$probability = DistributionValidations::validateProbability($probability);
$degrees = DistributionValidations::validateInt($degrees);
} catch (Exception $e) {
return $e->getMessage();
}
if ($degrees <= 0) {
return Functions::NAN();
}
$callback = function ($value) use ($degrees) {
return self::distribution($value, $degrees, 2);
};
$newtonRaphson = new NewtonRaphson($callback);
return $newtonRaphson->execute($probability);
}
/**
* @return float
*/
private static function calculateDistribution(float $value, int $degrees, int $tails)
{
// tdist, which finds the probability that corresponds to a given value
// of t with k degrees of freedom. This algorithm is translated from a
// pascal function on p81 of "Statistical Computing in Pascal" by D
// Cooke, A H Craven & G M Clark (1985: Edward Arnold (Pubs.) Ltd:
// London). The above Pascal algorithm is itself a translation of the
// fortran algoritm "AS 3" by B E Cooper of the Atlas Computer
// Laboratory as reported in (among other places) "Applied Statistics
// Algorithms", editied by P Griffiths and I D Hill (1985; Ellis
// Horwood Ltd.; W. Sussex, England).
$tterm = $degrees;
$ttheta = atan2($value, sqrt($tterm));
$tc = cos($ttheta);
$ts = sin($ttheta);
if (($degrees % 2) === 1) {
$ti = 3;
$tterm = $tc;
} else {
$ti = 2;
$tterm = 1;
}
$tsum = $tterm;
while ($ti < $degrees) {
$tterm *= $tc * $tc * ($ti - 1) / $ti;
$tsum += $tterm;
$ti += 2;
}
$tsum *= $ts;
if (($degrees % 2) == 1) {
$tsum = Functions::M_2DIVPI * ($tsum + $ttheta);
}
$tValue = 0.5 * (1 + $tsum);
if ($tails == 1) {
return 1 - abs($tValue);
}
return 1 - abs((1 - $tValue) - $tValue);
}
}
ChiSquared.php 0000644 00000022337 15053602601 0007312 0 ustar 00 getMessage();
}
if ($degrees < 1) {
return Functions::NAN();
}
if ($value < 0) {
if (Functions::getCompatibilityMode() == Functions::COMPATIBILITY_GNUMERIC) {
return 1;
}
return Functions::NAN();
}
return 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2) / Gamma::gammaValue($degrees / 2));
}
/**
* CHIDIST.
*
* Returns the one-tailed probability of the chi-squared distribution.
*
* @param mixed $value Float value for which we want the probability
* @param mixed $degrees Integer degrees of freedom
* @param mixed $cumulative Boolean value indicating if we want the cdf (true) or the pdf (false)
*
* @return float|string
*/
public static function distributionLeftTail($value, $degrees, $cumulative)
{
$value = Functions::flattenSingleValue($value);
$degrees = Functions::flattenSingleValue($degrees);
$cumulative = Functions::flattenSingleValue($cumulative);
try {
$value = DistributionValidations::validateFloat($value);
$degrees = DistributionValidations::validateInt($degrees);
$cumulative = DistributionValidations::validateBool($cumulative);
} catch (Exception $e) {
return $e->getMessage();
}
if ($degrees < 1) {
return Functions::NAN();
}
if ($value < 0) {
if (Functions::getCompatibilityMode() == Functions::COMPATIBILITY_GNUMERIC) {
return 1;
}
return Functions::NAN();
}
if ($cumulative === true) {
return 1 - self::distributionRightTail($value, $degrees);
}
return (($value ** (($degrees / 2) - 1) * exp(-$value / 2))) /
((2 ** ($degrees / 2)) * Gamma::gammaValue($degrees / 2));
}
/**
* CHIINV.
*
* Returns the inverse of the right-tailed probability of the chi-squared distribution.
*
* @param mixed $probability Float probability at which you want to evaluate the distribution
* @param mixed $degrees Integer degrees of freedom
*
* @return float|string
*/
public static function inverseRightTail($probability, $degrees)
{
$probability = Functions::flattenSingleValue($probability);
$degrees = Functions::flattenSingleValue($degrees);
try {
$probability = DistributionValidations::validateProbability($probability);
$degrees = DistributionValidations::validateInt($degrees);
} catch (Exception $e) {
return $e->getMessage();
}
if ($degrees < 1) {
return Functions::NAN();
}
$callback = function ($value) use ($degrees) {
return 1 - (Gamma::incompleteGamma($degrees / 2, $value / 2)
/ Gamma::gammaValue($degrees / 2));
};
$newtonRaphson = new NewtonRaphson($callback);
return $newtonRaphson->execute($probability);
}
/**
* CHIINV.
*
* Returns the inverse of the left-tailed probability of the chi-squared distribution.
*
* @param mixed $probability Float probability at which you want to evaluate the distribution
* @param mixed $degrees Integer degrees of freedom
*
* @return float|string
*/
public static function inverseLeftTail($probability, $degrees)
{
$probability = Functions::flattenSingleValue($probability);
$degrees = Functions::flattenSingleValue($degrees);
try {
$probability = DistributionValidations::validateProbability($probability);
$degrees = DistributionValidations::validateInt($degrees);
} catch (Exception $e) {
return $e->getMessage();
}
if ($degrees < 1) {
return Functions::NAN();
}
return self::inverseLeftTailCalculation($probability, $degrees);
}
/**
* CHITEST.
*
* Uses the chi-square test to calculate the probability that the differences between two supplied data sets
* (of observed and expected frequencies), are likely to be simply due to sampling error,
* or if they are likely to be real.
*
* @param mixed $actual an array of observed frequencies
* @param mixed $expected an array of expected frequencies
*
* @return float|string
*/
public static function test($actual, $expected)
{
$rows = count($actual);
$actual = Functions::flattenArray($actual);
$expected = Functions::flattenArray($expected);
$columns = count($actual) / $rows;
$countActuals = count($actual);
$countExpected = count($expected);
if ($countActuals !== $countExpected || $countActuals === 1) {
return Functions::NAN();
}
$result = 0.0;
for ($i = 0; $i < $countActuals; ++$i) {
if ($expected[$i] == 0.0) {
return Functions::DIV0();
} elseif ($expected[$i] < 0.0) {
return Functions::NAN();
}
$result += (($actual[$i] - $expected[$i]) ** 2) / $expected[$i];
}
$degrees = self::degrees($rows, $columns);
$result = self::distributionRightTail($result, $degrees);
return $result;
}
protected static function degrees(int $rows, int $columns): int
{
if ($rows === 1) {
return $columns - 1;
} elseif ($columns === 1) {
return $rows - 1;
}
return ($columns - 1) * ($rows - 1);
}
private static function inverseLeftTailCalculation(float $probability, int $degrees): float
{
// bracket the root
$min = 0;
$sd = sqrt(2.0 * $degrees);
$max = 2 * $sd;
$s = -1;
while ($s * self::pchisq($max, $degrees) > $probability * $s) {
$min = $max;
$max += 2 * $sd;
}
// Find root using bisection
$chi2 = 0.5 * ($min + $max);
while (($max - $min) > self::EPS * $chi2) {
if ($s * self::pchisq($chi2, $degrees) > $probability * $s) {
$min = $chi2;
} else {
$max = $chi2;
}
$chi2 = 0.5 * ($min + $max);
}
return $chi2;
}
private static function pchisq($chi2, $degrees)
{
return self::gammp($degrees, 0.5 * $chi2);
}
private static function gammp($n, $x)
{
if ($x < 0.5 * $n + 1) {
return self::gser($n, $x);
}
return 1 - self::gcf($n, $x);
}
// Return the incomplete gamma function P(n/2,x) evaluated by
// series representation. Algorithm from numerical recipe.
// Assume that n is a positive integer and x>0, won't check arguments.
// Relative error controlled by the eps parameter
private static function gser($n, $x)
{
$gln = Gamma::ln($n / 2);
$a = 0.5 * $n;
$ap = $a;
$sum = 1.0 / $a;
$del = $sum;
for ($i = 1; $i < 101; ++$i) {
++$ap;
$del = $del * $x / $ap;
$sum += $del;
if ($del < $sum * self::EPS) {
break;
}
}
return $sum * exp(-$x + $a * log($x) - $gln);
}
// Return the incomplete gamma function Q(n/2,x) evaluated by
// its continued fraction representation. Algorithm from numerical recipe.
// Assume that n is a postive integer and x>0, won't check arguments.
// Relative error controlled by the eps parameter
private static function gcf($n, $x)
{
$gln = Gamma::ln($n / 2);
$a = 0.5 * $n;
$b = $x + 1 - $a;
$fpmin = 1.e-300;
$c = 1 / $fpmin;
$d = 1 / $b;
$h = $d;
for ($i = 1; $i < 101; ++$i) {
$an = -$i * ($i - $a);
$b += 2;
$d = $an * $d + $b;
if (abs($d) < $fpmin) {
$d = $fpmin;
}
$c = $b + $an / $c;
if (abs($c) < $fpmin) {
$c = $fpmin;
}
$d = 1 / $d;
$del = $d * $c;
$h = $h * $del;
if (abs($del - 1) < self::EPS) {
break;
}
}
return $h * exp(-$x + $a * log($x) - $gln);
}
}
Fisher.php 0000644 00000003366 15053602601 0006503 0 ustar 00 getMessage();
}
if (($value <= -1) || ($value >= 1)) {
return Functions::NAN();
}
return 0.5 * log((1 + $value) / (1 - $value));
}
/**
* FISHERINV.
*
* Returns the inverse of the Fisher transformation. Use this transformation when
* analyzing correlations between ranges or arrays of data. If y = FISHER(x), then
* FISHERINV(y) = x.
*
* @param mixed $probability Float probability at which you want to evaluate the distribution
*
* @return float|string
*/
public static function inverse($probability)
{
$probability = Functions::flattenSingleValue($probability);
try {
DistributionValidations::validateFloat($probability);
} catch (Exception $e) {
return $e->getMessage();
}
return (exp(2 * $probability) - 1) / (exp(2 * $probability) + 1);
}
}
StandardNormal.php 0000644 00000007241 15053602601 0010170 0 ustar 00 getMessage();
}
if (($sampleSuccesses < 0) || ($sampleSuccesses > $sampleNumber) || ($sampleSuccesses > $populationSuccesses)) {
return Functions::NAN();
}
if (($sampleNumber <= 0) || ($sampleNumber > $populationNumber)) {
return Functions::NAN();
}
if (($populationSuccesses <= 0) || ($populationSuccesses > $populationNumber)) {
return Functions::NAN();
}
$successesPopulationAndSample = (float) Combinations::withoutRepetition($populationSuccesses, $sampleSuccesses);
$numbersPopulationAndSample = (float) Combinations::withoutRepetition($populationNumber, $sampleNumber);
$adjustedPopulationAndSample = (float) Combinations::withoutRepetition(
$populationNumber - $populationSuccesses,
$sampleNumber - $sampleSuccesses
);
return $successesPopulationAndSample * $adjustedPopulationAndSample / $numbersPopulationAndSample;
}
}
Weibull.php 0000644 00000003307 15053602601 0006661 0 ustar 00 getMessage();
}
if (($value < 0) || ($alpha <= 0) || ($beta <= 0)) {
return Functions::NAN();
}
if ($cumulative) {
return 1 - exp(0 - ($value / $beta) ** $alpha);
}
return ($alpha / $beta ** $alpha) * $value ** ($alpha - 1) * exp(0 - ($value / $beta) ** $alpha);
}
}
Beta.php 0000644 00000020615 15053602601 0006132 0 ustar 00 getMessage();
}
if ($rMin > $rMax) {
$tmp = $rMin;
$rMin = $rMax;
$rMax = $tmp;
}
if (($value < $rMin) || ($value > $rMax) || ($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax)) {
return Functions::NAN();
}
$value -= $rMin;
$value /= ($rMax - $rMin);
return self::incompleteBeta($value, $alpha, $beta);
}
/**
* BETAINV.
*
* Returns the inverse of the Beta distribution.
*
* @param mixed $probability Float probability at which you want to evaluate the distribution
* @param mixed $alpha Parameter to the distribution as a float
* @param mixed $beta Parameter to the distribution as a float
* @param mixed $rMin Minimum value as a float
* @param mixed $rMax Maximum value as a float
*
* @return float|string
*/
public static function inverse($probability, $alpha, $beta, $rMin = 0.0, $rMax = 1.0)
{
$probability = Functions::flattenSingleValue($probability);
$alpha = Functions::flattenSingleValue($alpha);
$beta = Functions::flattenSingleValue($beta);
$rMin = ($rMin === null) ? 0.0 : Functions::flattenSingleValue($rMin);
$rMax = ($rMax === null) ? 1.0 : Functions::flattenSingleValue($rMax);
try {
$probability = DistributionValidations::validateProbability($probability);
$alpha = DistributionValidations::validateFloat($alpha);
$beta = DistributionValidations::validateFloat($beta);
$rMax = DistributionValidations::validateFloat($rMax);
$rMin = DistributionValidations::validateFloat($rMin);
} catch (Exception $e) {
return $e->getMessage();
}
if ($rMin > $rMax) {
$tmp = $rMin;
$rMin = $rMax;
$rMax = $tmp;
}
if (($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax) || ($probability <= 0.0)) {
return Functions::NAN();
}
return self::calculateInverse($probability, $alpha, $beta, $rMin, $rMax);
}
/**
* @return float|string
*/
private static function calculateInverse(float $probability, float $alpha, float $beta, float $rMin, float $rMax)
{
$a = 0;
$b = 2;
$i = 0;
while ((($b - $a) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
$guess = ($a + $b) / 2;
$result = self::distribution($guess, $alpha, $beta);
if (($result === $probability) || ($result === 0.0)) {
$b = $a;
} elseif ($result > $probability) {
$b = $guess;
} else {
$a = $guess;
}
}
if ($i === self::MAX_ITERATIONS) {
return Functions::NA();
}
return round($rMin + $guess * ($rMax - $rMin), 12);
}
/**
* Incomplete beta function.
*
* @author Jaco van Kooten
* @author Paul Meagher
*
* The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992).
*
* @param float $x require 0<=x<=1
* @param float $p require p>0
* @param float $q require q>0
*
* @return float 0 if x<0, p<=0, q<=0 or p+q>2.55E305 and 1 if x>1 to avoid errors and over/underflow
*/
public static function incompleteBeta(float $x, float $p, float $q): float
{
if ($x <= 0.0) {
return 0.0;
} elseif ($x >= 1.0) {
return 1.0;
} elseif (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) {
return 0.0;
}
$beta_gam = exp((0 - self::logBeta($p, $q)) + $p * log($x) + $q * log(1.0 - $x));
if ($x < ($p + 1.0) / ($p + $q + 2.0)) {
return $beta_gam * self::betaFraction($x, $p, $q) / $p;
}
return 1.0 - ($beta_gam * self::betaFraction(1 - $x, $q, $p) / $q);
}
// Function cache for logBeta function
private static $logBetaCacheP = 0.0;
private static $logBetaCacheQ = 0.0;
private static $logBetaCacheResult = 0.0;
/**
* The natural logarithm of the beta function.
*
* @param float $p require p>0
* @param float $q require q>0
*
* @return float 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow
*
* @author Jaco van Kooten
*/
private static function logBeta(float $p, float $q): float
{
if ($p != self::$logBetaCacheP || $q != self::$logBetaCacheQ) {
self::$logBetaCacheP = $p;
self::$logBetaCacheQ = $q;
if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) {
self::$logBetaCacheResult = 0.0;
} else {
self::$logBetaCacheResult = Gamma::logGamma($p) + Gamma::logGamma($q) - Gamma::logGamma($p + $q);
}
}
return self::$logBetaCacheResult;
}
/**
* Evaluates of continued fraction part of incomplete beta function.
* Based on an idea from Numerical Recipes (W.H. Press et al, 1992).
*
* @author Jaco van Kooten
*/
private static function betaFraction(float $x, float $p, float $q): float
{
$c = 1.0;
$sum_pq = $p + $q;
$p_plus = $p + 1.0;
$p_minus = $p - 1.0;
$h = 1.0 - $sum_pq * $x / $p_plus;
if (abs($h) < self::XMININ) {
$h = self::XMININ;
}
$h = 1.0 / $h;
$frac = $h;
$m = 1;
$delta = 0.0;
while ($m <= self::MAX_ITERATIONS && abs($delta - 1.0) > Functions::PRECISION) {
$m2 = 2 * $m;
// even index for d
$d = $m * ($q - $m) * $x / (($p_minus + $m2) * ($p + $m2));
$h = 1.0 + $d * $h;
if (abs($h) < self::XMININ) {
$h = self::XMININ;
}
$h = 1.0 / $h;
$c = 1.0 + $d / $c;
if (abs($c) < self::XMININ) {
$c = self::XMININ;
}
$frac *= $h * $c;
// odd index for d
$d = -($p + $m) * ($sum_pq + $m) * $x / (($p + $m2) * ($p_plus + $m2));
$h = 1.0 + $d * $h;
if (abs($h) < self::XMININ) {
$h = self::XMININ;
}
$h = 1.0 / $h;
$c = 1.0 + $d / $c;
if (abs($c) < self::XMININ) {
$c = self::XMININ;
}
$delta = $h * $c;
$frac *= $delta;
++$m;
}
return $frac;
}
private static function betaValue(float $a, float $b): float
{
return (Gamma::gammaValue($a) * Gamma::gammaValue($b)) /
Gamma::gammaValue($a + $b);
}
private static function regularizedIncompleteBeta(float $value, float $a, float $b): float
{
return self::incompleteBeta($value, $a, $b) / self::betaValue($a, $b);
}
}
GammaBase.php 0000644 00000026213 15053602601 0007074 0 ustar 00 Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) {
// Apply Newton-Raphson step
$result = self::calculateDistribution($x, $alpha, $beta, true);
$error = $result - $probability;
if ($error == 0.0) {
$dx = 0;
} elseif ($error < 0.0) {
$xLo = $x;
} else {
$xHi = $x;
}
$pdf = self::calculateDistribution($x, $alpha, $beta, false);
// Avoid division by zero
if ($pdf !== 0.0) {
$dx = $error / $pdf;
$xNew = $x - $dx;
}
// If the NR fails to converge (which for example may be the
// case if the initial guess is too rough) we apply a bisection
// step to determine a more narrow interval around the root.
if (($xNew < $xLo) || ($xNew > $xHi) || ($pdf == 0.0)) {
$xNew = ($xLo + $xHi) / 2;
$dx = $xNew - $x;
}
$x = $xNew;
}
if ($i === self::MAX_ITERATIONS) {
return Functions::NA();
}
return $x;
}
//
// Implementation of the incomplete Gamma function
//
public static function incompleteGamma(float $a, float $x): float
{
static $max = 32;
$summer = 0;
for ($n = 0; $n <= $max; ++$n) {
$divisor = $a;
for ($i = 1; $i <= $n; ++$i) {
$divisor *= ($a + $i);
}
$summer += ($x ** $n / $divisor);
}
return $x ** $a * exp(0 - $x) * $summer;
}
//
// Implementation of the Gamma function
//
public static function gammaValue(float $value): float
{
if ($value == 0.0) {
return 0;
}
static $p0 = 1.000000000190015;
static $p = [
1 => 76.18009172947146,
2 => -86.50532032941677,
3 => 24.01409824083091,
4 => -1.231739572450155,
5 => 1.208650973866179e-3,
6 => -5.395239384953e-6,
];
$y = $x = $value;
$tmp = $x + 5.5;
$tmp -= ($x + 0.5) * log($tmp);
$summer = $p0;
for ($j = 1; $j <= 6; ++$j) {
$summer += ($p[$j] / ++$y);
}
return exp(0 - $tmp + log(self::SQRT2PI * $summer / $x));
}
/**
* logGamma function.
*
* @version 1.1
*
* @author Jaco van Kooten
*
* Original author was Jaco van Kooten. Ported to PHP by Paul Meagher.
*
* The natural logarithm of the gamma function.
* Based on public domain NETLIB (Fortran) code by W. J. Cody and L. Stoltz
* Applied Mathematics Division
* Argonne National Laboratory
* Argonne, IL 60439
*
* References:
*
* - W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for the Natural
* Logarithm of the Gamma Function,' Math. Comp. 21, 1967, pp. 198-203.
* - K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May, 1969.
* - Hart, Et. Al., Computer Approximations, Wiley and sons, New York, 1968.
*
*
*
* From the original documentation:
*
*
* This routine calculates the LOG(GAMMA) function for a positive real argument X.
* Computation is based on an algorithm outlined in references 1 and 2.
* The program uses rational functions that theoretically approximate LOG(GAMMA)
* to at least 18 significant decimal digits. The approximation for X > 12 is from
* reference 3, while approximations for X < 12.0 are similar to those in reference
* 1, but are unpublished. The accuracy achieved depends on the arithmetic system,
* the compiler, the intrinsic functions, and proper selection of the
* machine-dependent constants.
*
*
* Error returns:
* The program returns the value XINF for X .LE. 0.0 or when overflow would occur.
* The computation is believed to be free of underflow and overflow.
*
*
* @return float MAX_VALUE for x < 0.0 or when overflow would occur, i.e. x > 2.55E305
*/
// Log Gamma related constants
private const LG_D1 = -0.5772156649015328605195174;
private const LG_D2 = 0.4227843350984671393993777;
private const LG_D4 = 1.791759469228055000094023;
private const LG_P1 = [
4.945235359296727046734888,
201.8112620856775083915565,
2290.838373831346393026739,
11319.67205903380828685045,
28557.24635671635335736389,
38484.96228443793359990269,
26377.48787624195437963534,
7225.813979700288197698961,
];
private const LG_P2 = [
4.974607845568932035012064,
542.4138599891070494101986,
15506.93864978364947665077,
184793.2904445632425417223,
1088204.76946882876749847,
3338152.967987029735917223,
5106661.678927352456275255,
3074109.054850539556250927,
];
private const LG_P4 = [
14745.02166059939948905062,
2426813.369486704502836312,
121475557.4045093227939592,
2663432449.630976949898078,
29403789566.34553899906876,
170266573776.5398868392998,
492612579337.743088758812,
560625185622.3951465078242,
];
private const LG_Q1 = [
67.48212550303777196073036,
1113.332393857199323513008,
7738.757056935398733233834,
27639.87074403340708898585,
54993.10206226157329794414,
61611.22180066002127833352,
36351.27591501940507276287,
8785.536302431013170870835,
];
private const LG_Q2 = [
183.0328399370592604055942,
7765.049321445005871323047,
133190.3827966074194402448,
1136705.821321969608938755,
5267964.117437946917577538,
13467014.54311101692290052,
17827365.30353274213975932,
9533095.591844353613395747,
];
private const LG_Q4 = [
2690.530175870899333379843,
639388.5654300092398984238,
41355999.30241388052042842,
1120872109.61614794137657,
14886137286.78813811542398,
101680358627.2438228077304,
341747634550.7377132798597,
446315818741.9713286462081,
];
private const LG_C = [
-0.001910444077728,
8.4171387781295e-4,
-5.952379913043012e-4,
7.93650793500350248e-4,
-0.002777777777777681622553,
0.08333333333333333331554247,
0.0057083835261,
];
// Rough estimate of the fourth root of logGamma_xBig
private const LG_FRTBIG = 2.25e76;
private const PNT68 = 0.6796875;
// Function cache for logGamma
private static $logGammaCacheResult = 0.0;
private static $logGammaCacheX = 0.0;
public static function logGamma(float $x): float
{
if ($x == self::$logGammaCacheX) {
return self::$logGammaCacheResult;
}
$y = $x;
if ($y > 0.0 && $y <= self::LOG_GAMMA_X_MAX_VALUE) {
if ($y <= self::EPS) {
$res = -log($y);
} elseif ($y <= 1.5) {
$res = self::logGamma1($y);
} elseif ($y <= 4.0) {
$res = self::logGamma2($y);
} elseif ($y <= 12.0) {
$res = self::logGamma3($y);
} else {
$res = self::logGamma4($y);
}
} else {
// --------------------------
// Return for bad arguments
// --------------------------
$res = self::MAX_VALUE;
}
// ------------------------------
// Final adjustments and return
// ------------------------------
self::$logGammaCacheX = $x;
self::$logGammaCacheResult = $res;
return $res;
}
private static function logGamma1(float $y)
{
// ---------------------
// EPS .LT. X .LE. 1.5
// ---------------------
if ($y < self::PNT68) {
$corr = -log($y);
$xm1 = $y;
} else {
$corr = 0.0;
$xm1 = $y - 1.0;
}
$xden = 1.0;
$xnum = 0.0;
if ($y <= 0.5 || $y >= self::PNT68) {
for ($i = 0; $i < 8; ++$i) {
$xnum = $xnum * $xm1 + self::LG_P1[$i];
$xden = $xden * $xm1 + self::LG_Q1[$i];
}
return $corr + $xm1 * (self::LG_D1 + $xm1 * ($xnum / $xden));
}
$xm2 = $y - 1.0;
for ($i = 0; $i < 8; ++$i) {
$xnum = $xnum * $xm2 + self::LG_P2[$i];
$xden = $xden * $xm2 + self::LG_Q2[$i];
}
return $corr + $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
}
private static function logGamma2(float $y)
{
// ---------------------
// 1.5 .LT. X .LE. 4.0
// ---------------------
$xm2 = $y - 2.0;
$xden = 1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; ++$i) {
$xnum = $xnum * $xm2 + self::LG_P2[$i];
$xden = $xden * $xm2 + self::LG_Q2[$i];
}
return $xm2 * (self::LG_D2 + $xm2 * ($xnum / $xden));
}
protected static function logGamma3(float $y)
{
// ----------------------
// 4.0 .LT. X .LE. 12.0
// ----------------------
$xm4 = $y - 4.0;
$xden = -1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; ++$i) {
$xnum = $xnum * $xm4 + self::LG_P4[$i];
$xden = $xden * $xm4 + self::LG_Q4[$i];
}
return self::LG_D4 + $xm4 * ($xnum / $xden);
}
protected static function logGamma4(float $y)
{
// ---------------------------------
// Evaluate for argument .GE. 12.0
// ---------------------------------
$res = 0.0;
if ($y <= self::LG_FRTBIG) {
$res = self::LG_C[6];
$ysq = $y * $y;
for ($i = 0; $i < 6; ++$i) {
$res = $res / $ysq + self::LG_C[$i];
}
$res /= $y;
$corr = log($y);
$res = $res + log(self::SQRT2PI) - 0.5 * $corr;
$res += $y * ($corr - 1.0);
}
return $res;
}
}
Exponential.php 0000644 00000003146 15053602601 0007545 0 ustar 00 getMessage();
}
if (($value < 0) || ($lambda < 0)) {
return Functions::NAN();
}
if ($cumulative === true) {
return 1 - exp(0 - $value * $lambda);
}
return $lambda * exp(0 - $value * $lambda);
}
}
LogNormal.php 0000644 00000010357 15053602601 0007153 0 ustar 00 getMessage();
}
if (($value <= 0) || ($stdDev <= 0)) {
return Functions::NAN();
}
return StandardNormal::cumulative((log($value) - $mean) / $stdDev);
}
/**
* LOGNORM.DIST.
*
* Returns the lognormal distribution of x, where ln(x) is normally distributed
* with parameters mean and standard_dev.
*
* @param mixed $value Float value for which we want the probability
* @param mixed $mean Mean value as a float
* @param mixed $stdDev Standard Deviation as a float
* @param mixed $cumulative Boolean value indicating if we want the cdf (true) or the pdf (false)
*
* @return float|string The result, or a string containing an error
*/
public static function distribution($value, $mean, $stdDev, $cumulative = false)
{
$value = Functions::flattenSingleValue($value);
$mean = Functions::flattenSingleValue($mean);
$stdDev = Functions::flattenSingleValue($stdDev);
$cumulative = Functions::flattenSingleValue($cumulative);
try {
$value = DistributionValidations::validateFloat($value);
$mean = DistributionValidations::validateFloat($mean);
$stdDev = DistributionValidations::validateFloat($stdDev);
$cumulative = DistributionValidations::validateBool($cumulative);
} catch (Exception $e) {
return $e->getMessage();
}
if (($value <= 0) || ($stdDev <= 0)) {
return Functions::NAN();
}
if ($cumulative === true) {
return StandardNormal::distribution((log($value) - $mean) / $stdDev, true);
}
return (1 / (sqrt(2 * M_PI) * $stdDev * $value)) *
exp(0 - ((log($value) - $mean) ** 2 / (2 * $stdDev ** 2)));
}
/**
* LOGINV.
*
* Returns the inverse of the lognormal cumulative distribution
*
* @param mixed $probability Float probability for which we want the value
* @param mixed $mean Mean Value as a float
* @param mixed $stdDev Standard Deviation as a float
*
* @return float|string The result, or a string containing an error
*
* @TODO Try implementing P J Acklam's refinement algorithm for greater
* accuracy if I can get my head round the mathematics
* (as described at) http://home.online.no/~pjacklam/notes/invnorm/
*/
public static function inverse($probability, $mean, $stdDev)
{
$probability = Functions::flattenSingleValue($probability);
$mean = Functions::flattenSingleValue($mean);
$stdDev = Functions::flattenSingleValue($stdDev);
try {
$probability = DistributionValidations::validateProbability($probability);
$mean = DistributionValidations::validateFloat($mean);
$stdDev = DistributionValidations::validateFloat($stdDev);
} catch (Exception $e) {
return $e->getMessage();
}
if ($stdDev <= 0) {
return Functions::NAN();
}
return exp($mean + $stdDev * StandardNormal::inverse($probability));
}
}
Gamma.php 0000644 00000007365 15053602601 0006310 0 ustar 00 getMessage();
}
if ((((int) $value) == ((float) $value)) && $value <= 0.0) {
return Functions::NAN();
}
return self::gammaValue($value);
}
/**
* GAMMADIST.
*
* Returns the gamma distribution.
*
* @param mixed $value Float Value at which you want to evaluate the distribution
* @param mixed $a Parameter to the distribution as a float
* @param mixed $b Parameter to the distribution as a float
* @param mixed $cumulative Boolean value indicating if we want the cdf (true) or the pdf (false)
*
* @return float|string
*/
public static function distribution($value, $a, $b, $cumulative)
{
$value = Functions::flattenSingleValue($value);
$a = Functions::flattenSingleValue($a);
$b = Functions::flattenSingleValue($b);
try {
$value = DistributionValidations::validateFloat($value);
$a = DistributionValidations::validateFloat($a);
$b = DistributionValidations::validateFloat($b);
$cumulative = DistributionValidations::validateBool($cumulative);
} catch (Exception $e) {
return $e->getMessage();
}
if (($value < 0) || ($a <= 0) || ($b <= 0)) {
return Functions::NAN();
}
return self::calculateDistribution($value, $a, $b, $cumulative);
}
/**
* GAMMAINV.
*
* Returns the inverse of the Gamma distribution.
*
* @param mixed $probability Float probability at which you want to evaluate the distribution
* @param mixed $alpha Parameter to the distribution as a float
* @param mixed $beta Parameter to the distribution as a float
*
* @return float|string
*/
public static function inverse($probability, $alpha, $beta)
{
$probability = Functions::flattenSingleValue($probability);
$alpha = Functions::flattenSingleValue($alpha);
$beta = Functions::flattenSingleValue($beta);
try {
$probability = DistributionValidations::validateProbability($probability);
$alpha = DistributionValidations::validateFloat($alpha);
$beta = DistributionValidations::validateFloat($beta);
} catch (Exception $e) {
return $e->getMessage();
}
if (($alpha <= 0.0) || ($beta <= 0.0)) {
return Functions::NAN();
}
return self::calculateInverse($probability, $alpha, $beta);
}
/**
* GAMMALN.
*
* Returns the natural logarithm of the gamma function.
*
* @param mixed $value Float Value at which you want to evaluate the distribution
*
* @return float|string
*/
public static function ln($value)
{
$value = Functions::flattenSingleValue($value);
try {
$value = DistributionValidations::validateFloat($value);
} catch (Exception $e) {
return $e->getMessage();
}
if ($value <= 0) {
return Functions::NAN();
}
return log(self::gammaValue($value));
}
}
Binomial.php 0000644 00000017544 15053602601 0007020 0 ustar 00 getMessage();
}
if (($value < 0) || ($value > $trials)) {
return Functions::NAN();
}
if ($cumulative) {
return self::calculateCumulativeBinomial($value, $trials, $probability);
}
return Combinations::withoutRepetition($trials, $value) * $probability ** $value
* (1 - $probability) ** ($trials - $value);
}
/**
* BINOM.DIST.RANGE.
*
* Returns returns the Binomial Distribution probability for the number of successes from a specified number
* of trials falling into a specified range.
*
* @param mixed $trials Integer number of trials
* @param mixed $probability Probability of success on each trial as a float
* @param mixed $successes The integer number of successes in trials
* @param mixed $limit Upper limit for successes in trials as null, or an integer
* If null, then this will indicate the same as the number of Successes
*
* @return float|string
*/
public static function range($trials, $probability, $successes, $limit = null)
{
$trials = Functions::flattenSingleValue($trials);
$probability = Functions::flattenSingleValue($probability);
$successes = Functions::flattenSingleValue($successes);
$limit = ($limit === null) ? $successes : Functions::flattenSingleValue($limit);
try {
$trials = DistributionValidations::validateInt($trials);
$probability = DistributionValidations::validateProbability($probability);
$successes = DistributionValidations::validateInt($successes);
$limit = DistributionValidations::validateInt($limit);
} catch (Exception $e) {
return $e->getMessage();
}
if (($successes < 0) || ($successes > $trials)) {
return Functions::NAN();
}
if (($limit < 0) || ($limit > $trials) || $limit < $successes) {
return Functions::NAN();
}
$summer = 0;
for ($i = $successes; $i <= $limit; ++$i) {
$summer += Combinations::withoutRepetition($trials, $i) * $probability ** $i
* (1 - $probability) ** ($trials - $i);
}
return $summer;
}
/**
* NEGBINOMDIST.
*
* Returns the negative binomial distribution. NEGBINOMDIST returns the probability that
* there will be number_f failures before the number_s-th success, when the constant
* probability of a success is probability_s. This function is similar to the binomial
* distribution, except that the number of successes is fixed, and the number of trials is
* variable. Like the binomial, trials are assumed to be independent.
*
* @param mixed $failures Number of Failures as an integer
* @param mixed $successes Threshold number of Successes as an integer
* @param mixed $probability Probability of success on each trial as a float
*
* @return float|string The result, or a string containing an error
*
* TODO Add support for the cumulative flag not present for NEGBINOMDIST, but introduced for NEGBINOM.DIST
* The cumulative default should be false to reflect the behaviour of NEGBINOMDIST
*/
public static function negative($failures, $successes, $probability)
{
$failures = Functions::flattenSingleValue($failures);
$successes = Functions::flattenSingleValue($successes);
$probability = Functions::flattenSingleValue($probability);
try {
$failures = DistributionValidations::validateInt($failures);
$successes = DistributionValidations::validateInt($successes);
$probability = DistributionValidations::validateProbability($probability);
} catch (Exception $e) {
return $e->getMessage();
}
if (($failures < 0) || ($successes < 1)) {
return Functions::NAN();
}
if (Functions::getCompatibilityMode() == Functions::COMPATIBILITY_GNUMERIC) {
if (($failures + $successes - 1) <= 0) {
return Functions::NAN();
}
}
return (Combinations::withoutRepetition($failures + $successes - 1, $successes - 1))
* ($probability ** $successes) * ((1 - $probability) ** $failures);
}
/**
* CRITBINOM.
*
* Returns the smallest value for which the cumulative binomial distribution is greater
* than or equal to a criterion value
*
* @param mixed $trials number of Bernoulli trials as an integer
* @param mixed $probability probability of a success on each trial as a float
* @param mixed $alpha criterion value as a float
*
* @return int|string
*/
public static function inverse($trials, $probability, $alpha)
{
$trials = Functions::flattenSingleValue($trials);
$probability = Functions::flattenSingleValue($probability);
$alpha = Functions::flattenSingleValue($alpha);
try {
$trials = DistributionValidations::validateInt($trials);
$probability = DistributionValidations::validateProbability($probability);
$alpha = DistributionValidations::validateFloat($alpha);
} catch (Exception $e) {
return $e->getMessage();
}
if ($trials < 0) {
return Functions::NAN();
} elseif (($alpha < 0.0) || ($alpha > 1.0)) {
return Functions::NAN();
}
$successes = 0;
while ($successes <= $trials) {
$result = self::calculateCumulativeBinomial($successes, $trials, $probability);
if ($result >= $alpha) {
break;
}
++$successes;
}
return $successes;
}
/**
* @return float|int
*/
private static function calculateCumulativeBinomial(int $value, int $trials, float $probability)
{
$summer = 0;
for ($i = 0; $i <= $value; ++$i) {
$summer += Combinations::withoutRepetition($trials, $i) * $probability ** $i
* (1 - $probability) ** ($trials - $i);
}
return $summer;
}
}
F.php 0000644 00000004132 15053602601 0005440 0 ustar 00 getMessage();
}
if ($value < 0 || $u < 1 || $v < 1) {
return Functions::NAN();
}
if ($cumulative) {
$adjustedValue = ($u * $value) / ($u * $value + $v);
return Beta::incompleteBeta($adjustedValue, $u / 2, $v / 2);
}
return (Gamma::gammaValue(($v + $u) / 2) /
(Gamma::gammaValue($u / 2) * Gamma::gammaValue($v / 2))) *
(($u / $v) ** ($u / 2)) *
(($value ** (($u - 2) / 2)) / ((1 + ($u / $v) * $value) ** (($u + $v) / 2)));
}
}
DistributionValidations.php 0000644 00000001230 15053602601 0012124 0 ustar 00 1.0) {
throw new Exception(Functions::NAN());
}
return $probability;
}
}
Normal.php 0000644 00000014175 15053602601 0006513 0 ustar 00 getMessage();
}
if ($stdDev < 0) {
return Functions::NAN();
}
if ($cumulative) {
return 0.5 * (1 + Engineering\Erf::erfValue(($value - $mean) / ($stdDev * sqrt(2))));
}
return (1 / (self::SQRT2PI * $stdDev)) * exp(0 - (($value - $mean) ** 2 / (2 * ($stdDev * $stdDev))));
}
/**
* NORMINV.
*
* Returns the inverse of the normal cumulative distribution for the specified mean and standard deviation.
*
* @param mixed $probability Float probability for which we want the value
* @param mixed $mean Mean Value as a float
* @param mixed $stdDev Standard Deviation as a float
*
* @return float|string The result, or a string containing an error
*/
public static function inverse($probability, $mean, $stdDev)
{
$probability = Functions::flattenSingleValue($probability);
$mean = Functions::flattenSingleValue($mean);
$stdDev = Functions::flattenSingleValue($stdDev);
try {
$probability = DistributionValidations::validateProbability($probability);
$mean = DistributionValidations::validateFloat($mean);
$stdDev = DistributionValidations::validateFloat($stdDev);
} catch (Exception $e) {
return $e->getMessage();
}
if ($stdDev < 0) {
return Functions::NAN();
}
return (self::inverseNcdf($probability) * $stdDev) + $mean;
}
/*
* inverse_ncdf.php
* -------------------
* begin : Friday, January 16, 2004
* copyright : (C) 2004 Michael Nickerson
* email : nickersonm@yahoo.com
*
*/
private static function inverseNcdf($p)
{
// Inverse ncdf approximation by Peter J. Acklam, implementation adapted to
// PHP by Michael Nickerson, using Dr. Thomas Ziegler's C implementation as
// a guide. http://home.online.no/~pjacklam/notes/invnorm/index.html
// I have not checked the accuracy of this implementation. Be aware that PHP
// will truncate the coeficcients to 14 digits.
// You have permission to use and distribute this function freely for
// whatever purpose you want, but please show common courtesy and give credit
// where credit is due.
// Input paramater is $p - probability - where 0 < p < 1.
// Coefficients in rational approximations
static $a = [
1 => -3.969683028665376e+01,
2 => 2.209460984245205e+02,
3 => -2.759285104469687e+02,
4 => 1.383577518672690e+02,
5 => -3.066479806614716e+01,
6 => 2.506628277459239e+00,
];
static $b = [
1 => -5.447609879822406e+01,
2 => 1.615858368580409e+02,
3 => -1.556989798598866e+02,
4 => 6.680131188771972e+01,
5 => -1.328068155288572e+01,
];
static $c = [
1 => -7.784894002430293e-03,
2 => -3.223964580411365e-01,
3 => -2.400758277161838e+00,
4 => -2.549732539343734e+00,
5 => 4.374664141464968e+00,
6 => 2.938163982698783e+00,
];
static $d = [
1 => 7.784695709041462e-03,
2 => 3.224671290700398e-01,
3 => 2.445134137142996e+00,
4 => 3.754408661907416e+00,
];
// Define lower and upper region break-points.
$p_low = 0.02425; //Use lower region approx. below this
$p_high = 1 - $p_low; //Use upper region approx. above this
if (0 < $p && $p < $p_low) {
// Rational approximation for lower region.
$q = sqrt(-2 * log($p));
return ((((($c[1] * $q + $c[2]) * $q + $c[3]) * $q + $c[4]) * $q + $c[5]) * $q + $c[6]) /
(((($d[1] * $q + $d[2]) * $q + $d[3]) * $q + $d[4]) * $q + 1);
} elseif ($p_high < $p && $p < 1) {
// Rational approximation for upper region.
$q = sqrt(-2 * log(1 - $p));
return -((((($c[1] * $q + $c[2]) * $q + $c[3]) * $q + $c[4]) * $q + $c[5]) * $q + $c[6]) /
(((($d[1] * $q + $d[2]) * $q + $d[3]) * $q + $d[4]) * $q + 1);
}
// Rational approximation for central region.
$q = $p - 0.5;
$r = $q * $q;
return ((((($a[1] * $r + $a[2]) * $r + $a[3]) * $r + $a[4]) * $r + $a[5]) * $r + $a[6]) * $q /
((((($b[1] * $r + $b[2]) * $r + $b[3]) * $r + $b[4]) * $r + $b[5]) * $r + 1);
}
}
NewtonRaphson.php 0000644 00000003067 15053602601 0010066 0 ustar 00 callback = $callback;
}
public function execute(float $probability)
{
$xLo = 100;
$xHi = 0;
$dx = 1;
$x = $xNew = 1;
$i = 0;
while ((abs($dx) > Functions::PRECISION) && ($i++ < self::MAX_ITERATIONS)) {
// Apply Newton-Raphson step
$result = call_user_func($this->callback, $x);
$error = $result - $probability;
if ($error == 0.0) {
$dx = 0;
} elseif ($error < 0.0) {
$xLo = $x;
} else {
$xHi = $x;
}
// Avoid division by zero
if ($result != 0.0) {
$dx = $error / $result;
$xNew = $x - $dx;
}
// If the NR fails to converge (which for example may be the
// case if the initial guess is too rough) we apply a bisection
// step to determine a more narrow interval around the root.
if (($xNew < $xLo) || ($xNew > $xHi) || ($result == 0.0)) {
$xNew = ($xLo + $xHi) / 2;
$dx = $xNew - $x;
}
$x = $xNew;
}
if ($i == self::MAX_ITERATIONS) {
return Functions::NA();
}
return $x;
}
}
Poisson.php 0000644 00000003370 15053602601 0006710 0 ustar 00 getMessage();
}
if (($value < 0) || ($mean < 0)) {
return Functions::NAN();
}
if ($cumulative) {
$summer = 0;
$floor = floor($value);
for ($i = 0; $i <= $floor; ++$i) {
$summer += $mean ** $i / MathTrig\Factorial::fact($i);
}
return exp(0 - $mean) * $summer;
}
return (exp(0 - $mean) * $mean ** $value) / MathTrig\Factorial::fact($value);
}
}