Skip to content

Commit 47fa8ae

Browse files
gburtinigburtini
authored
Initial pass of cleanup. (#29)
Linting, minor bug fixes, refactoring in small places. Mostly nitpicky and PHP-Code_Sniffer stuff for PSR2 standard.
1 parent 12cfa24 commit 47fa8ae

25 files changed

+2723
-2627
lines changed

‎README.md‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -42,8 +42,6 @@ All supported distributions are in the namespace `gburtini\Distributions` and im
4242
`gburtini\Distributions` contains the distribution classes as indicated above.
4343
`gburtini\Distributions\Accessories` contains BetaFunction and GammaFunction, two classes containing accessory functions for computing complete, incomplete and inverse beta and gamma functions numerically.
4444

45-
While unsupported, if you are using a version of PHP pre-namespaces, the `ugly/` directory implements pseudonamespaced ("ugly") names, with the prefix GBPDP\_ and can be used directly. If you have access to namespaces (PHP 5.3+) you should use the composer-compatible namespacing to interact with the classes. You **really** should upgrade your PHP version to a [supported version](http://php.net/supported-versions.php).
46-
4745
## Example
4846

4947
Examples are provided in a comment at the top of most of the implementation files. In general, you should be able to use the parametrization listed above under "Supported Distributions" to create classes that implement the methods under "Interfaces".

‎src/gburtini/Distributions/Accessories/BetaFunction.php‎

Lines changed: 150 additions & 149 deletions
Original file line numberDiff line numberDiff line change
@@ -1,150 +1,151 @@
11
<?php
2-
namespace gburtini\Distributions\Accessories;
3-
4-
require_once dirname(__FILE__) . "/GammaFunction.php";
5-
6-
// these accesory functions are "necessary" to compute many beta distribution statistics, but we don't actually use them (here, they are used in the T distribution calculations currently) because the Gamma shortcut calculator is faster.
7-
class BetaFunction {
8-
public static function betaFunction($x, $y) {
9-
return exp(GammaFunction::logGammaFunction($x) + GammaFunction::logGammaFunction($y) - GammaFunction::logGammaFunction($x+$y));
10-
}
11-
12-
public static function inverseIncompleteBetaFunction($p, $a, $b) {
13-
// from jStat.ibetainv
14-
$EPS = 1e-8;
15-
$a1 = $a - 1;
16-
$b1 = $b - 1;
17-
$j = 0;
18-
//$lna; $lnb; $pp; $t; $u; $err; $x; $al; $h; $w; $afac;
19-
if ($p <= 0)
20-
return 0;
21-
if ($p >= 1)
22-
return 1;
23-
if ($a >= 1 && $b >= 1) {
24-
$pp = ($p < 0.5) ? $p : 1 - $p;
25-
$t = sqrt(-2 * log($pp));
26-
$x = (2.30753 + $t * 0.27061) / (1 + $t* (0.99229 + $t * 0.04481)) - $t;
27-
if ($p < 0.5)
28-
$x = -$x;
29-
$al = ($x * $x - 3) / 6;
30-
$h = 2 / (1 / (2 * $a - 1) + 1 / (2 * $b - 1));
31-
$w = ($x * sqrt($al + $h) / $h) - (1 / (2 * $b - 1) - 1 / (2 * $a - 1)) * ($al + 5 / 6 - 2 / (3 * $h));
32-
$x = $a / ($a + $b * exp(2 * $w));
33-
} else {
34-
$lna = log($a / ($a + $b));
35-
$lnb = log($b / ($a + $b));
36-
$t = exp($a * $lna) / $a;
37-
$u = exp($b * $lnb) / $b;
38-
$w = $t + $u;
39-
if ($p < $t / $w)
40-
$x = pow($a * $w * $p, 1 / $a);
41-
else
42-
$x = 1 - pow($b * $w * (1 - $p), 1 / $b);
43-
}
44-
$afac = -GammaFunction::logGammaFunction($a) - GammaFunction::logGammaFunction($b) + GammaFunction::logGammaFunction($a + $b);
45-
for(; $j < 10; $j++) {
46-
if ($x === 0 || $x === 1)
47-
return $x;
48-
$err = self::incompleteBetaFunction($x, $a, $b) - $p;
49-
$t = exp($a1 * log($x) + $b1 * log(1 - $x) + $afac);
50-
$u = $err / $t;
51-
$x -= ($t = $u / (1 - 0.5 * min(1, $u * ($a1 / $x - $b1 / (1 - $x)))));
52-
if ($x <= 0)
53-
$x = 0.5 * ($x + $t);
54-
if ($x >= 1)
55-
$x = 0.5 * ($x + $t + 1);
56-
if (abs($t) < $EPS * $x && $j > 0)
57-
break;
58-
}
59-
return $x;
60-
}
61-
62-
public static function incompleteBetaFunction($x, $a, $b) {
63-
$g_ab = GammaFunction::logGammaFunction($a+$b);
64-
$g_a = GammaFunction::logGammaFunction($a);
65-
$g_b = GammaFunction::logGammaFunction($b);
66-
67-
$bt = exp($g_ab - $g_a - $g_b + $a*log($x)+$b*log(1.0-$x));
68-
69-
if ($x == 0)
70-
$bt = 0;
71-
72-
if ($x < (($a + 1.0)/($a + $b + 2.0)))
73-
return $bt* self::continuedFraction($a,$b,$x)/$a;
74-
else
75-
return 1 - ($bt*self::continuedFraction($b,$a,1.0-$x)/$b);
76-
}
77-
78-
// see http://functions.wolfram.com/GammaBetaErf/Beta3/10/
79-
public static function continuedFraction($a, $b, $x) {
80-
$maxit = 100;
81-
$eps = 3e-16;
82-
$fpmin = 1e-30;
83-
/*
84-
// mentioning these causes E_NOTICEs in HipHop.
85-
$aa;
86-
$c;
87-
$d;
88-
$del;
89-
$h;
90-
$qab;
91-
$qam;
92-
$qap;
93-
*/
94-
95-
$qab = $a + $b;
96-
$qap = $a + 1;
97-
$qam = $a - 1;
98-
99-
$c = 1.0;
100-
$d = 1.0 - $qab*$x/$qap;
101-
102-
if (abs($d)<$fpmin)
103-
$d = $fpmin;
104-
105-
$d = 1.0/$d;
106-
107-
$h = $d;
108-
109-
//$m2;
110-
111-
for ($m = 1; $m < $maxit; $m++)
112-
{
113-
$m2 = 2*$m;
114-
$aa = $m*($b-$m)*$x/(($qam + $m2)*($a + $m2));
115-
$d = 1.0 + $aa*$d;
116-
117-
if (abs($d)<$fpmin)
118-
$d = $fpmin;
119-
120-
$c = 1.0 + $aa/$c;
121-
122-
if (abs($c)<$fpmin)
123-
$c = $fpmin;
124-
125-
$d = 1.0/$d;
126-
$h = $h*$d*$c;
127-
$aa = -($a + $m)*($qab + $m)*$x/(($a+$m2)*($qap+$m2));
128-
$d = 1.0 + $aa*$d;
129-
130-
if (abs($d)<$fpmin)
131-
$d = $fpmin;
132-
133-
$c = 1.0 + $aa/$c;
134-
135-
if (abs($c)<$fpmin)
136-
$c = $fpmin;
137-
138-
$d = 1.0/$d;
139-
$del = $d*$c;
140-
$h = $h*$del;
141-
142-
if (abs($del-1.0)< $eps)
143-
{
144-
break;
145-
}
146-
}
147-
return $h;
148-
}
149-
150-
}
2+
namespace gburtini\Distributions\Accessories;
3+
4+
class BetaFunction
5+
{
6+
public static function betaFunction($x, $y)
7+
{
8+
return exp(GammaFunction::logGammaFunction($x) + GammaFunction::logGammaFunction($y) - GammaFunction::logGammaFunction($x + $y));
9+
}
10+
11+
public static function inverseIncompleteBetaFunction($p, $a, $b)
12+
{
13+
// from jStat.ibetainv
14+
$EPS = 1e-8;
15+
$a1 = $a - 1;
16+
$b1 = $b - 1;
17+
$j = 0;
18+
//$lna; $lnb; $pp; $t; $u; $err; $x; $al; $h; $w; $afac;
19+
if ($p <= 0) {
20+
return 0;
21+
}
22+
if ($p >= 1) {
23+
return 1;
24+
}
25+
if ($a >= 1 && $b >= 1) {
26+
$pp = ($p < 0.5) ? $p : 1 - $p;
27+
$t = sqrt(-2 * log($pp));
28+
$x = (2.30753 + $t * 0.27061) / (1 + $t* (0.99229 + $t * 0.04481)) - $t;
29+
if ($p < 0.5) {
30+
$x = -$x;
31+
}
32+
$al = ($x * $x - 3) / 6;
33+
$h = 2 / (1 / (2 * $a - 1) + 1 / (2 * $b - 1));
34+
$w = ($x * sqrt($al + $h) / $h) - (1 / (2 * $b - 1) - 1 / (2 * $a - 1)) * ($al + 5 / 6 - 2 / (3 * $h));
35+
$x = $a / ($a + $b * exp(2 * $w));
36+
} else {
37+
$lna = log($a / ($a + $b));
38+
$lnb = log($b / ($a + $b));
39+
$t = exp($a * $lna) / $a;
40+
$u = exp($b * $lnb) / $b;
41+
$w = $t + $u;
42+
if ($p < $t / $w) {
43+
$x = pow($a * $w * $p, 1 / $a);
44+
} else {
45+
$x = 1 - pow($b * $w * (1 - $p), 1 / $b);
46+
}
47+
}
48+
$afac = -GammaFunction::logGammaFunction($a) - GammaFunction::logGammaFunction($b) + GammaFunction::logGammaFunction($a + $b);
49+
for (; $j < 10; $j++) {
50+
if ($x === 0 || $x === 1) {
51+
return $x;
52+
}
53+
$err = self::incompleteBetaFunction($x, $a, $b) - $p;
54+
$t = exp($a1 * log($x) + $b1 * log(1 - $x) + $afac);
55+
$u = $err / $t;
56+
$x -= ($t = $u / (1 - 0.5 * min(1, $u * ($a1 / $x - $b1 / (1 - $x)))));
57+
if ($x <= 0) {
58+
$x = 0.5 * ($x + $t);
59+
}
60+
if ($x >= 1) {
61+
$x = 0.5 * ($x + $t + 1);
62+
}
63+
if (abs($t) < $EPS * $x && $j > 0) {
64+
break;
65+
}
66+
}
67+
return $x;
68+
}
69+
70+
public static function incompleteBetaFunction($x, $a, $b)
71+
{
72+
$g_ab = GammaFunction::logGammaFunction($a+$b);
73+
$g_a = GammaFunction::logGammaFunction($a);
74+
$g_b = GammaFunction::logGammaFunction($b);
75+
76+
$bt = exp($g_ab - $g_a - $g_b + $a*log($x)+$b*log(1.0-$x));
77+
78+
if ($x == 0) {
79+
$bt = 0;
80+
}
81+
82+
if ($x < (($a + 1.0)/($a + $b + 2.0))) {
83+
return $bt* self::continuedFraction($a, $b, $x)/$a;
84+
} else {
85+
return 1 - ($bt*self::continuedFraction($b, $a, 1.0-$x)/$b);
86+
}
87+
}
88+
89+
// see http://functions.wolfram.com/GammaBetaErf/Beta3/10/
90+
public static function continuedFraction($a, $b, $x)
91+
{
92+
$maxit = 100;
93+
$eps = 3e-16;
94+
$fpmin = 1e-30;
95+
96+
$qab = $a + $b;
97+
$qap = $a + 1;
98+
$qam = $a - 1;
99+
100+
$c = 1.0;
101+
$d = 1.0 - $qab*$x/$qap;
102+
103+
if (abs($d)<$fpmin) {
104+
$d = $fpmin;
105+
}
106+
107+
$d = 1.0/$d;
108+
109+
$h = $d;
110+
111+
for ($m = 1; $m < $maxit; $m++) {
112+
$m2 = 2*$m;
113+
$aa = $m*($b-$m)*$x/(($qam + $m2)*($a + $m2));
114+
$d = 1.0 + $aa*$d;
115+
116+
if (abs($d)<$fpmin) {
117+
$d = $fpmin;
118+
}
119+
120+
$c = 1.0 + $aa/$c;
121+
122+
if (abs($c)<$fpmin) {
123+
$c = $fpmin;
124+
}
125+
126+
$d = 1.0/$d;
127+
$h = $h*$d*$c;
128+
$aa = -($a + $m)*($qab + $m)*$x/(($a+$m2)*($qap+$m2));
129+
$d = 1.0 + $aa*$d;
130+
131+
if (abs($d)<$fpmin) {
132+
$d = $fpmin;
133+
}
134+
135+
$c = 1.0 + $aa/$c;
136+
137+
if (abs($c)<$fpmin) {
138+
$c = $fpmin;
139+
}
140+
141+
$d = 1.0/$d;
142+
$del = $d*$c;
143+
$h = $h*$del;
144+
145+
if (abs($del-1.0)< $eps) {
146+
break;
147+
}
148+
}
149+
return $h;
150+
}
151+
}

0 commit comments

Comments
 (0)