tablepress
/
libraries
/
vendor
/
PhpSpreadsheet
/
Calculation
/
Statistical
/
Distributions
/
Beta.php
Beta.php in TablePress – Tables in WordPress made easy 3.4, at libraries/vendor/PhpSpreadsheet/Calculation/Statistical/Distributions/Beta.php
| 1 | <?php |
| 2 | |
| 3 | namespace TablePress\PhpOffice\PhpSpreadsheet\Calculation\Statistical\Distributions; |
| 4 | |
| 5 | use TablePress\PhpOffice\PhpSpreadsheet\Calculation\ArrayEnabled; |
| 6 | use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Exception; |
| 7 | use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Functions; |
| 8 | use TablePress\PhpOffice\PhpSpreadsheet\Calculation\Information\ExcelError; |
| 9 | |
| 10 | class Beta |
| 11 | { |
| 12 | use ArrayEnabled; |
| 13 | |
| 14 | private const MAX_ITERATIONS = 256; |
| 15 | |
| 16 | private const LOG_GAMMA_X_MAX_VALUE = 2.55e305; |
| 17 | |
| 18 | private const XMININ = 2.23e-308; |
| 19 | |
| 20 | /** |
| 21 | * BETADIST. |
| 22 | * |
| 23 | * Returns the beta distribution. |
| 24 | * |
| 25 | * @param mixed $value Float value at which you want to evaluate the distribution |
| 26 | * Or can be an array of values |
| 27 | * @param mixed $alpha Parameter to the distribution as a float |
| 28 | * Or can be an array of values |
| 29 | * @param mixed $beta Parameter to the distribution as a float |
| 30 | * Or can be an array of values |
| 31 | * @param mixed $rMin as a float |
| 32 | * Or can be an array of values |
| 33 | * @param mixed $rMax as a float |
| 34 | * Or can be an array of values |
| 35 | * |
| 36 | * @return array<mixed>|float|string If an array of numbers is passed as an argument, then the returned result will also be an array |
| 37 | * with the same dimensions |
| 38 | */ |
| 39 | public static function distribution($value, $alpha, $beta, $rMin = 0.0, $rMax = 1.0) |
| 40 | { |
| 41 | if (is_array($value) || is_array($alpha) || is_array($beta) || is_array($rMin) || is_array($rMax)) { |
| 42 | return self::evaluateArrayArguments([self::class, __FUNCTION__], $value, $alpha, $beta, $rMin, $rMax); |
| 43 | } |
| 44 | |
| 45 | $rMin = $rMin ?? 0.0; |
| 46 | $rMax = $rMax ?? 1.0; |
| 47 | |
| 48 | try { |
| 49 | $value = DistributionValidations::validateFloat($value); |
| 50 | $alpha = DistributionValidations::validateFloat($alpha); |
| 51 | $beta = DistributionValidations::validateFloat($beta); |
| 52 | $rMax = DistributionValidations::validateFloat($rMax); |
| 53 | $rMin = DistributionValidations::validateFloat($rMin); |
| 54 | } catch (Exception $e) { |
| 55 | return $e->getMessage(); |
| 56 | } |
| 57 | |
| 58 | if ($rMin > $rMax) { |
| 59 | $tmp = $rMin; |
| 60 | $rMin = $rMax; |
| 61 | $rMax = $tmp; |
| 62 | } |
| 63 | if (($value < $rMin) || ($value > $rMax) || ($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax)) { |
| 64 | return ExcelError::NAN(); |
| 65 | } |
| 66 | |
| 67 | $value -= $rMin; |
| 68 | $value /= ($rMax - $rMin); |
| 69 | |
| 70 | return self::incompleteBeta($value, $alpha, $beta); |
| 71 | } |
| 72 | |
| 73 | /** |
| 74 | * BETAINV. |
| 75 | * |
| 76 | * Returns the inverse of the Beta distribution. |
| 77 | * |
| 78 | * @param mixed $probability Float probability at which you want to evaluate the distribution |
| 79 | * Or can be an array of values |
| 80 | * @param mixed $alpha Parameter to the distribution as a float |
| 81 | * Or can be an array of values |
| 82 | * @param mixed $beta Parameter to the distribution as a float |
| 83 | * Or can be an array of values |
| 84 | * @param mixed $rMin Minimum value as a float |
| 85 | * Or can be an array of values |
| 86 | * @param mixed $rMax Maximum value as a float |
| 87 | * Or can be an array of values |
| 88 | * |
| 89 | * @return array<mixed>|float|string If an array of numbers is passed as an argument, then the returned result will also be an array |
| 90 | * with the same dimensions |
| 91 | */ |
| 92 | public static function inverse($probability, $alpha, $beta, $rMin = 0.0, $rMax = 1.0) |
| 93 | { |
| 94 | if (is_array($probability) || is_array($alpha) || is_array($beta) || is_array($rMin) || is_array($rMax)) { |
| 95 | return self::evaluateArrayArguments([self::class, __FUNCTION__], $probability, $alpha, $beta, $rMin, $rMax); |
| 96 | } |
| 97 | |
| 98 | $rMin = $rMin ?? 0.0; |
| 99 | $rMax = $rMax ?? 1.0; |
| 100 | |
| 101 | try { |
| 102 | $probability = DistributionValidations::validateProbability($probability); |
| 103 | $alpha = DistributionValidations::validateFloat($alpha); |
| 104 | $beta = DistributionValidations::validateFloat($beta); |
| 105 | $rMax = DistributionValidations::validateFloat($rMax); |
| 106 | $rMin = DistributionValidations::validateFloat($rMin); |
| 107 | } catch (Exception $e) { |
| 108 | return $e->getMessage(); |
| 109 | } |
| 110 | |
| 111 | if ($rMin > $rMax) { |
| 112 | $tmp = $rMin; |
| 113 | $rMin = $rMax; |
| 114 | $rMax = $tmp; |
| 115 | } |
| 116 | if (($alpha <= 0) || ($beta <= 0) || ($rMin == $rMax) || ($probability <= 0.0)) { |
| 117 | return ExcelError::NAN(); |
| 118 | } |
| 119 | if (($alpha + $beta) > self::LOG_GAMMA_X_MAX_VALUE) { |
| 120 | // incompleteBeta declines to evaluate here and returns 0 for every x, |
| 121 | // so there is no quantile to search for. |
| 122 | return ExcelError::NAN(); |
| 123 | } |
| 124 | |
| 125 | return self::calculateInverse($probability, $alpha, $beta, $rMin, $rMax); |
| 126 | } |
| 127 | |
| 128 | /** |
| 129 | * @return string|float |
| 130 | */ |
| 131 | private static function calculateInverse(float $probability, float $alpha, float $beta, float $rMin, float $rMax) |
| 132 | { |
| 133 | $a = 0; |
| 134 | $b = 2; |
| 135 | $guess = ($a + $b) / 2; |
| 136 | |
| 137 | $i = 0; |
| 138 | while ((($b - $a) > Functions::PRECISION) && (++$i <= self::MAX_ITERATIONS)) { |
| 139 | $guess = ($a + $b) / 2; |
| 140 | $result = self::distribution($guess, $alpha, $beta); |
| 141 | if ($result === $probability) { |
| 142 | $b = $a; |
| 143 | } elseif ($result > $probability) { |
| 144 | $b = $guess; |
| 145 | } else { |
| 146 | $a = $guess; |
| 147 | } |
| 148 | } |
| 149 | |
| 150 | if ($i === self::MAX_ITERATIONS) { |
| 151 | return ExcelError::NA(); |
| 152 | } |
| 153 | |
| 154 | return round($rMin + $guess * ($rMax - $rMin), 12); |
| 155 | } |
| 156 | |
| 157 | /** |
| 158 | * Incomplete beta function. |
| 159 | * |
| 160 | * @author Jaco van Kooten |
| 161 | * @author Paul Meagher |
| 162 | * |
| 163 | * The computation is based on formulas from Numerical Recipes, Chapter 6.4 (W.H. Press et al, 1992). |
| 164 | * |
| 165 | * @param float $x require 0<=x<=1 |
| 166 | * @param float $p require p>0 |
| 167 | * @param float $q require q>0 |
| 168 | * |
| 169 | * @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 |
| 170 | */ |
| 171 | public static function incompleteBeta(float $x, float $p, float $q): float |
| 172 | { |
| 173 | if ($x <= 0.0) { |
| 174 | return 0.0; |
| 175 | } elseif ($x >= 1.0) { |
| 176 | return 1.0; |
| 177 | } elseif (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) { |
| 178 | return 0.0; |
| 179 | } |
| 180 | |
| 181 | $beta_gam = exp((0 - self::logBeta($p, $q)) + $p * log($x) + $q * log(1.0 - $x)); |
| 182 | if ($x < ($p + 1.0) / ($p + $q + 2.0)) { |
| 183 | return $beta_gam * self::betaFraction($x, $p, $q) / $p; |
| 184 | } |
| 185 | |
| 186 | return 1.0 - ($beta_gam * self::betaFraction(1 - $x, $q, $p) / $q); |
| 187 | } |
| 188 | |
| 189 | // Function cache for logBeta function |
| 190 | |
| 191 | private static float $logBetaCacheP = 0.0; |
| 192 | |
| 193 | private static float $logBetaCacheQ = 0.0; |
| 194 | |
| 195 | private static float $logBetaCacheResult = 0.0; |
| 196 | |
| 197 | /** |
| 198 | * The natural logarithm of the beta function. |
| 199 | * |
| 200 | * @param float $p require p>0 |
| 201 | * @param float $q require q>0 |
| 202 | * |
| 203 | * @return float 0 if p<=0, q<=0 or p+q>2.55E305 to avoid errors and over/underflow |
| 204 | * |
| 205 | * @author Jaco van Kooten |
| 206 | */ |
| 207 | private static function logBeta(float $p, float $q): float |
| 208 | { |
| 209 | if ($p != self::$logBetaCacheP || $q != self::$logBetaCacheQ) { |
| 210 | self::$logBetaCacheP = $p; |
| 211 | self::$logBetaCacheQ = $q; |
| 212 | if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > self::LOG_GAMMA_X_MAX_VALUE)) { |
| 213 | self::$logBetaCacheResult = 0.0; |
| 214 | } else { |
| 215 | self::$logBetaCacheResult = Gamma::logGamma($p) + Gamma::logGamma($q) - Gamma::logGamma($p + $q); |
| 216 | } |
| 217 | } |
| 218 | |
| 219 | return self::$logBetaCacheResult; |
| 220 | } |
| 221 | |
| 222 | /** |
| 223 | * Evaluates of continued fraction part of incomplete beta function. |
| 224 | * Based on an idea from Numerical Recipes (W.H. Press et al, 1992). |
| 225 | * |
| 226 | * @author Jaco van Kooten |
| 227 | */ |
| 228 | private static function betaFraction(float $x, float $p, float $q): float |
| 229 | { |
| 230 | $c = 1.0; |
| 231 | $sum_pq = $p + $q; |
| 232 | $p_plus = $p + 1.0; |
| 233 | $p_minus = $p - 1.0; |
| 234 | $h = 1.0 - $sum_pq * $x / $p_plus; |
| 235 | if (abs($h) < self::XMININ) { |
| 236 | $h = self::XMININ; |
| 237 | } |
| 238 | $h = 1.0 / $h; |
| 239 | $frac = $h; |
| 240 | $m = 1; |
| 241 | $delta = 0.0; |
| 242 | while ($m <= self::MAX_ITERATIONS && abs($delta - 1.0) > Functions::PRECISION) { |
| 243 | $m2 = 2 * $m; |
| 244 | // even index for d |
| 245 | $d = $m * ($q - $m) * $x / (($p_minus + $m2) * ($p + $m2)); |
| 246 | $h = 1.0 + $d * $h; |
| 247 | if (abs($h) < self::XMININ) { |
| 248 | $h = self::XMININ; |
| 249 | } |
| 250 | $h = 1.0 / $h; |
| 251 | $c = 1.0 + $d / $c; |
| 252 | if (abs($c) < self::XMININ) { |
| 253 | $c = self::XMININ; |
| 254 | } |
| 255 | $frac *= $h * $c; |
| 256 | // odd index for d |
| 257 | $d = -($p + $m) * ($sum_pq + $m) * $x / (($p + $m2) * ($p_plus + $m2)); |
| 258 | $h = 1.0 + $d * $h; |
| 259 | if (abs($h) < self::XMININ) { |
| 260 | $h = self::XMININ; |
| 261 | } |
| 262 | $h = 1.0 / $h; |
| 263 | $c = 1.0 + $d / $c; |
| 264 | if (abs($c) < self::XMININ) { |
| 265 | $c = self::XMININ; |
| 266 | } |
| 267 | $delta = $h * $c; |
| 268 | $frac *= $delta; |
| 269 | ++$m; |
| 270 | } |
| 271 | |
| 272 | return $frac; |
| 273 | } |
| 274 | |
| 275 | /* |
| 276 | private static function betaValue(float $a, float $b): float |
| 277 | { |
| 278 | return (Gamma::gammaValue($a) * Gamma::gammaValue($b)) / |
| 279 | Gamma::gammaValue($a + $b); |
| 280 | } |
| 281 | |
| 282 | private static function regularizedIncompleteBeta(float $value, float $a, float $b): float |
| 283 | { |
| 284 | return self::incompleteBeta($value, $a, $b) / self::betaValue($a, $b); |
| 285 | } |
| 286 | */ |
| 287 | } |
| 288 |