Skip to content

Instantly share code, notes, and snippets.

@kurokikaze
Created October 6, 2010 10:39
Show Gist options
  • Save kurokikaze/613154 to your computer and use it in GitHub Desktop.
Save kurokikaze/613154 to your computer and use it in GitHub Desktop.
<?php
/**
* Taken from PHPExcel
*
* Copyright (c) 2006 - 2009 PHPExcel
*
* This library is free software; you can redistribute it and/or
* modify it under the terms of the GNU Lesser General Public
* License as published by the Free Software Foundation; either
* version 2.1 of the License, or (at your option) any later version.
*
* This library 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
* Lesser General Public License for more details.
*
* You should have received a copy of the GNU Lesser General Public
* License along with this library; if not, write to the Free Software
* Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
*
* @category PHPExcel
* @package PHPExcel_Calculation
* @copyright Copyright (c) 2006 - 2009 PHPExcel (http://www.codeplex.com/PHPExcel)
* @license http://www.gnu.org/licenses/old-licenses/lgpl-2.1.txt LGPL
* @version 1.6.7, 2009-04-22
*/
define('EPS', 2.22e-16);
define('MAX_VALUE', 1.2e308);
define('LOG_GAMMA_X_MAX_VALUE', 2.55e305);
define('SQRT2PI', 2.5066282746310005024157652848110452530069867406099);
define('XMININ', 2.23e-308);
define('MAX_ITERATIONS', 150);
define('PRECISION', 8.88E-8);
define('EULER', 2.71828182845904523536);
class calculator {
private static $_errorCodes = array( 'null' => '#NULL!',
'divisionbyzero' => '#DIV/0!',
'value' => '#VALUE!',
'reference' => '#REF!',
'name' => '#NAME?',
'num' => '#NUM!',
'na' => '#N/A',
'gettingdata' => '#GETTING_DATA'
);
public static function flattenSingleValue($value = '') {
if (is_array($value)) {
$value = self::flattenSingleValue(array_pop($value));
}
return $value;
}
public static function BETAINV($probability, $alpha, $beta, $rMin=0, $rMax=1) {
$probability = self::flattenSingleValue($probability);
$alpha = self::flattenSingleValue($alpha);
$beta = self::flattenSingleValue($beta);
$rMin = self::flattenSingleValue($rMin);
$rMax = self::flattenSingleValue($rMax);
if ((is_numeric($probability)) && (is_numeric($alpha)) && (is_numeric($beta)) && (is_numeric($rMin)) && (is_numeric($rMax))) {
if (($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax) || ($probability <= 0) || ($probability > 1)) {
return self::$_errorCodes['num'];
}
if ($rMin > $rMax) {
$tmp = $rMin;
$rMin = $rMax;
$rMax = $tmp;
}
$a = 0;
$b = 2;
$maxIteration = 100;
$i = 0;
while ((($b - $a) > PRECISION) && ($i++ < MAX_ITERATIONS)) {
$guess = ($a + $b) / 2;
$result = self::BETADIST($guess, $alpha, $beta);
if (($result == $probability) || ($result == 0)) {
$b = $a;
} elseif ($result > $probability) {
$b = $guess;
} else {
$a = $guess;
}
}
if ($i == MAX_ITERATIONS) {
return self::$_errorCodes['na'];
}
return round($rMin + $guess * ($rMax - $rMin),12);
}
return self::$_errorCodes['value'];
}
public static function BETADIST($value,$alpha,$beta,$rMin=0,$rMax=1) {
$value = self::flattenSingleValue($value);
$alpha = self::flattenSingleValue($alpha);
$beta = self::flattenSingleValue($beta);
$rMin = self::flattenSingleValue($rMin);
$rMax = self::flattenSingleValue($rMax);
if ((is_numeric($value)) && (is_numeric($alpha)) && (is_numeric($beta)) && (is_numeric($rMin)) && (is_numeric($rMax))) {
if (($value < $rMin) || ($value > $rMax) || ($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax)) {
return self::$_errorCodes['num'];
}
if ($rMin > $rMax) {
$tmp = $rMin;
$rMin = $rMax;
$rMax = $tmp;
}
$value -= $rMin;
$value /= ($rMax - $rMin);
return self::incompleteBeta($value,$alpha,$beta);
}
return self::$_errorCodes['value'];
}
private static function betaFraction($x, $p, $q) {
$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) < XMININ) {
$h = XMININ;
}
$h = 1.0 / $h;
$frac = $h;
$m = 1;
$delta = 0.0;
while ($m <= MAX_ITERATIONS && abs($delta-1.0) > 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) < XMININ) {
$h = XMININ;
}
$h = 1.0 / $h;
$c = 1.0 + $d / $c;
if (abs($c) < XMININ) {
$c = 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) < XMININ) {
$h = XMININ;
}
$h = 1.0 / $h;
$c = 1.0 + $d / $c;
if (abs($c) < XMININ) {
$c = XMININ;
}
$delta = $h * $c;
$frac *= $delta;
++$m;
}
return $frac;
}
private static function incompleteBeta($x, $p, $q) {
if ($x <= 0.0) {
return 0.0;
} elseif ($x >= 1.0) {
return 1.0;
} elseif (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > 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;
} else {
return 1.0 - ($beta_gam * self::betaFraction(1 - $x, $q, $p) / $q);
}
}
private static $logBetaCache_p = 0.0;
private static $logBetaCache_q = 0.0;
private static $logBetaCache_result = 0.0;
private static function logBeta($p, $q) {
if ($p != self::$logBetaCache_p || $q != self::$logBetaCache_q) {
self::$logBetaCache_p = $p;
self::$logBetaCache_q = $q;
if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > LOG_GAMMA_X_MAX_VALUE)) {
self::$logBetaCache_result = 0.0;
} else {
self::$logBetaCache_result = self::logGamma($p) + self::logGamma($q) - self::logGamma($p + $q);
}
}
return self::$logBetaCache_result;
}
private static $logGammaCache_result = 0.0;
private static $logGammaCache_x = 0.0;
private static function logGamma($x) {
// Log Gamma related constants
static $lg_d1 = -0.5772156649015328605195174;
static $lg_d2 = 0.4227843350984671393993777;
static $lg_d4 = 1.791759469228055000094023;
static $lg_p1 = array( 4.945235359296727046734888,
201.8112620856775083915565,
2290.838373831346393026739,
11319.67205903380828685045,
28557.24635671635335736389,
38484.96228443793359990269,
26377.48787624195437963534,
7225.813979700288197698961 );
static $lg_p2 = array( 4.974607845568932035012064,
542.4138599891070494101986,
15506.93864978364947665077,
184793.2904445632425417223,
1088204.76946882876749847,
3338152.967987029735917223,
5106661.678927352456275255,
3074109.054850539556250927 );
static $lg_p4 = array( 14745.02166059939948905062,
2426813.369486704502836312,
121475557.4045093227939592,
2663432449.630976949898078,
29403789566.34553899906876,
170266573776.5398868392998,
492612579337.743088758812,
560625185622.3951465078242 );
static $lg_q1 = array( 67.48212550303777196073036,
1113.332393857199323513008,
7738.757056935398733233834,
27639.87074403340708898585,
54993.10206226157329794414,
61611.22180066002127833352,
36351.27591501940507276287,
8785.536302431013170870835 );
static $lg_q2 = array( 183.0328399370592604055942,
7765.049321445005871323047,
133190.3827966074194402448,
1136705.821321969608938755,
5267964.117437946917577538,
13467014.54311101692290052,
17827365.30353274213975932,
9533095.591844353613395747 );
static $lg_q4 = array( 2690.530175870899333379843,
639388.5654300092398984238,
41355999.30241388052042842,
1120872109.61614794137657,
14886137286.78813811542398,
101680358627.2438228077304,
341747634550.7377132798597,
446315818741.9713286462081 );
static $lg_c = array( -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
static $lg_frtbig = 2.25e76;
static $pnt68 = 0.6796875;
if ($x == self::$logGammaCache_x) {
return self::$logGammaCache_result;
}
$y = $x;
if ($y > 0.0 && $y <= LOG_GAMMA_X_MAX_VALUE) {
if ($y <= EPS) {
$res = -log(y);
} elseif ($y <= 1.5) {
// ---------------------
// EPS .LT. X .LE. 1.5
// ---------------------
if ($y < $pnt68) {
$corr = -log($y);
$xm1 = $y;
} else {
$corr = 0.0;
$xm1 = $y - 1.0;
}
if ($y <= 0.5 || $y >= $pnt68) {
$xden = 1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; ++$i) {
$xnum = $xnum * $xm1 + $lg_p1[$i];
$xden = $xden * $xm1 + $lg_q1[$i];
}
$res = $corr + $xm1 * ($lg_d1 + $xm1 * ($xnum / $xden));
} else {
$xm2 = $y - 1.0;
$xden = 1.0;
$xnum = 0.0;
for ($i = 0; $i < 8; ++$i) {
$xnum = $xnum * $xm2 + $lg_p2[$i];
$xden = $xden * $xm2 + $lg_q2[$i];
}
$res = $corr + $xm2 * ($lg_d2 + $xm2 * ($xnum / $xden));
}
} elseif ($y <= 4.0) {
// ---------------------
// 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 + $lg_p2[$i];
$xden = $xden * $xm2 + $lg_q2[$i];
}
$res = $xm2 * ($lg_d2 + $xm2 * ($xnum / $xden));
} elseif ($y <= 12.0) {
// ----------------------
// 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 + $lg_p4[$i];
$xden = $xden * $xm4 + $lg_q4[$i];
}
$res = $lg_d4 + $xm4 * ($xnum / $xden);
} else {
// ---------------------------------
// Evaluate for argument .GE. 12.0
// ---------------------------------
$res = 0.0;
if ($y <= $lg_frtbig) {
$res = $lg_c[6];
$ysq = $y * $y;
for ($i = 0; $i < 6; ++$i)
$res = $res / $ysq + $lg_c[$i];
}
$res /= $y;
$corr = log($y);
$res = $res + log(SQRT2PI) - 0.5 * $corr;
$res += $y * ($corr - 1.0);
}
} else {
// --------------------------
// Return for bad arguments
// --------------------------
$res = MAX_VALUE;
}
// ------------------------------
// Final adjustments and return
// ------------------------------
self::$logGammaCache_x = $x;
self::$logGammaCache_result = $res;
return $res;
}
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment