summaryrefslogtreecommitdiff
path: root/src/de/lmu/ifi/dbs/elki/math/statistics/distribution/GammaDistribution.java
diff options
context:
space:
mode:
Diffstat (limited to 'src/de/lmu/ifi/dbs/elki/math/statistics/distribution/GammaDistribution.java')
-rw-r--r--src/de/lmu/ifi/dbs/elki/math/statistics/distribution/GammaDistribution.java317
1 files changed, 132 insertions, 185 deletions
diff --git a/src/de/lmu/ifi/dbs/elki/math/statistics/distribution/GammaDistribution.java b/src/de/lmu/ifi/dbs/elki/math/statistics/distribution/GammaDistribution.java
index 21eebc51..1b9e2b42 100644
--- a/src/de/lmu/ifi/dbs/elki/math/statistics/distribution/GammaDistribution.java
+++ b/src/de/lmu/ifi/dbs/elki/math/statistics/distribution/GammaDistribution.java
@@ -4,7 +4,7 @@ package de.lmu.ifi.dbs.elki.math.statistics.distribution;
This file is part of ELKI:
Environment for Developing KDD-Applications Supported by Index-Structures
- Copyright (C) 2012
+ Copyright (C) 2013
Ludwig-Maximilians-Universität München
Lehr- und Forschungseinheit für Datenbanksysteme
ELKI Development Team
@@ -34,7 +34,7 @@ import de.lmu.ifi.dbs.elki.utilities.documentation.Reference;
*
* @author Erich Schubert
*/
-public class GammaDistribution implements DistributionWithRandom {
+public class GammaDistribution implements Distribution {
/**
* Euler–Mascheroni constant
*/
@@ -51,11 +51,22 @@ public class GammaDistribution implements DistributionWithRandom {
static final double[] LANCZOS = { 0.99999999999999709182, 57.156235665862923517, -59.597960355475491248, 14.136097974741747174, -0.49191381609762019978, .33994649984811888699e-4, .46523628927048575665e-4, -.98374475304879564677e-4, .15808870322491248884e-3, -.21026444172410488319e-3, .21743961811521264320e-3, -.16431810653676389022e-3, .84418223983852743293e-4, -.26190838401581408670e-4, .36899182659531622704e-5, };
/**
- * Numerical precision to use
+ * Numerical precision to use (data type dependent!)
+ *
+ * If you change this, make sure to test exhaustively!
*/
static final double NUM_PRECISION = 1E-15;
/**
+ * Maximum number of iterations for regularizedGammaP. To prevent degeneration
+ * for extreme values.
+ *
+ * FIXME: is this too high, too low? Can we improve behavior for extreme
+ * cases?
+ */
+ static final int MAX_ITERATIONS = 1000;
+
+ /**
* Alpha == k
*/
private final double k;
@@ -79,8 +90,8 @@ public class GammaDistribution implements DistributionWithRandom {
*/
public GammaDistribution(double k, double theta, Random random) {
super();
- if(k <= 0.0 || theta <= 0.0) {
- throw new IllegalArgumentException("Invalid parameters for Gamma distribution.");
+ if (!(k > 0.0) || !(theta > 0.0)) { // Note: also tests for NaNs!
+ throw new IllegalArgumentException("Invalid parameters for Gamma distribution: " + k + " " + theta);
}
this.k = k;
@@ -151,6 +162,9 @@ public class GammaDistribution implements DistributionWithRandom {
* @return cdf value
*/
public static double cdf(double val, double k, double theta) {
+ if (val < 0) {
+ return 0.;
+ }
return regularizedGammaP(k, val * theta);
}
@@ -163,6 +177,9 @@ public class GammaDistribution implements DistributionWithRandom {
* @return cdf value
*/
public static double logcdf(double val, double k, double theta) {
+ if (val < 0) {
+ return Double.NEGATIVE_INFINITY;
+ }
return logregularizedGammaP(k, val * theta);
}
@@ -175,18 +192,17 @@ public class GammaDistribution implements DistributionWithRandom {
* @return probability density
*/
public static double pdf(double x, double k, double theta) {
- if(x < 0) {
+ if (x < 0) {
return 0.0;
}
- if(x == 0) {
- if(k == 1.0) {
+ if (x == 0) {
+ if (k == 1.0) {
return theta;
- }
- else {
+ } else {
return 0.0;
}
}
- if(k == 1.0) {
+ if (k == 1.0) {
return Math.exp(-x * theta) * theta;
}
@@ -202,18 +218,17 @@ public class GammaDistribution implements DistributionWithRandom {
* @return probability density
*/
public static double logpdf(double x, double k, double theta) {
- if(x < 0) {
+ if (x < 0) {
return Double.NEGATIVE_INFINITY;
}
- if(x == 0) {
- if(k == 1.0) {
+ if (x == 0) {
+ if (k == 1.0) {
return Math.log(theta);
- }
- else {
+ } else {
return Double.NEGATIVE_INFINITY;
}
}
- if(k == 1.0) {
+ if (k == 1.0) {
return Math.log(theta) - x * theta;
}
@@ -232,14 +247,14 @@ public class GammaDistribution implements DistributionWithRandom {
* @return log(&#915;(x))
*/
public static double logGamma(final double x) {
- if(Double.isNaN(x) || (x <= 0.0)) {
+ if (Double.isNaN(x) || (x <= 0.0)) {
return Double.NaN;
}
double g = 607.0 / 128.0;
double tmp = x + g + .5;
tmp = (x + 0.5) * Math.log(tmp) - tmp;
double ser = LANCZOS[0];
- for(int i = LANCZOS.length - 1; i > 0; --i) {
+ for (int i = LANCZOS.length - 1; i > 0; --i) {
ser += LANCZOS[i] / (x + i);
}
return tmp + Math.log(MathUtil.SQRTTWOPI * ser / x);
@@ -275,30 +290,30 @@ public class GammaDistribution implements DistributionWithRandom {
*/
public static double regularizedGammaP(final double a, final double x) {
// Special cases
- if(Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
+ if (Double.isInfinite(a) || Double.isInfinite(x) || !(a > 0.0) || !(x >= 0.0)) {
return Double.NaN;
}
- if(x == 0.0) {
+ if (x == 0.0) {
return 0.0;
}
- if(x >= a + 1) {
+ if (x >= a + 1) {
// Expected to converge faster
return 1.0 - regularizedGammaQ(a, x);
}
// Loosely following "Numerical Recipes"
- double del = 1.0 / a;
- double sum = del;
- for(int n = 1; n < Integer.MAX_VALUE; n++) {
+ double term = 1.0 / a;
+ double sum = term;
+ for (int n = 1; n < MAX_ITERATIONS; n++) {
// compute next element in the series
- del *= x / (a + n);
- sum = sum + del;
- if(Math.abs(del / sum) < NUM_PRECISION || sum >= Double.POSITIVE_INFINITY) {
+ term = x / (a + n) * term;
+ sum = sum + term;
+ if (sum == Double.POSITIVE_INFINITY) {
+ return 1.0;
+ }
+ if (Math.abs(term / sum) < NUM_PRECISION) {
break;
}
}
- if(Double.isInfinite(sum)) {
- return 1.0;
- }
return Math.exp(-x + (a * Math.log(x)) - logGamma(a)) * sum;
}
@@ -316,13 +331,13 @@ public class GammaDistribution implements DistributionWithRandom {
*/
public static double logregularizedGammaP(final double a, final double x) {
// Special cases
- if(Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
+ if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
return Double.NaN;
}
- if(x == 0.0) {
+ if (x == 0.0) {
return Double.NEGATIVE_INFINITY;
}
- if(x >= a + 1) {
+ if (x >= a + 1) {
// Expected to converge faster
// FIXME: and in log?
return Math.log(1.0 - regularizedGammaQ(a, x));
@@ -330,15 +345,15 @@ public class GammaDistribution implements DistributionWithRandom {
// Loosely following "Numerical Recipes"
double del = 1.0 / a;
double sum = del;
- for(int n = 1; n < Integer.MAX_VALUE; n++) {
+ for (int n = 1; n < Integer.MAX_VALUE; n++) {
// compute next element in the series
del *= x / (a + n);
sum = sum + del;
- if(Math.abs(del / sum) < NUM_PRECISION || sum >= Double.POSITIVE_INFINITY) {
+ if (Math.abs(del / sum) < NUM_PRECISION || sum >= Double.POSITIVE_INFINITY) {
break;
}
}
- if(Double.isInfinite(sum)) {
+ if (Double.isInfinite(sum)) {
return 0;
}
// TODO: reread numerical recipes, can we replace log(sum)?
@@ -360,13 +375,13 @@ public class GammaDistribution implements DistributionWithRandom {
* @return Result
*/
public static double regularizedGammaQ(final double a, final double x) {
- if(Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
+ if (Double.isNaN(a) || Double.isNaN(x) || (a <= 0.0) || (x < 0.0)) {
return Double.NaN;
}
- if(x == 0.0) {
+ if (x == 0.0) {
return 1.0;
}
- if(x < a + 1.0) {
+ if (x < a + 1.0) {
// Expected to converge faster
return 1.0 - regularizedGammaP(a, x);
}
@@ -376,21 +391,21 @@ public class GammaDistribution implements DistributionWithRandom {
double c = 1.0 / FPMIN;
double d = 1.0 / b;
double fac = d;
- for(int i = 1; i < Integer.MAX_VALUE; i++) {
+ for (int i = 1; i < MAX_ITERATIONS; i++) {
double an = i * (a - i);
b += 2;
d = an * d + b;
- if(Math.abs(d) < FPMIN) {
+ if (Math.abs(d) < FPMIN) {
d = FPMIN;
}
c = b + an / c;
- if(Math.abs(c) < FPMIN) {
+ if (Math.abs(c) < FPMIN) {
c = FPMIN;
}
d = 1 / d;
double del = d * c;
fac *= del;
- if(Math.abs(del - 1.0) <= NUM_PRECISION) {
+ if (Math.abs(del - 1.0) <= NUM_PRECISION) {
break;
}
}
@@ -424,33 +439,30 @@ public class GammaDistribution implements DistributionWithRandom {
final double e4 = 0.041664508, e5 = 0.008345522, e6 = 0.001353826;
final double e7 = 0.000247453;
- if(k < 1.0) { // Base case, for small k
+ if (k < 1.0) { // Base case, for small k
final double b = 1.0 + 0.36788794412 * k; // Step 1
- while(true) {
+ while (true) {
final double p = b * random.nextDouble();
- if(p <= 1.0) { // when gds <= 1
+ if (p <= 1.0) { // when gds <= 1
final double gds = Math.exp(Math.log(p) / k);
- if(Math.log(random.nextDouble()) <= -gds) {
+ if (Math.log(random.nextDouble()) <= -gds) {
return (gds / theta);
}
- }
- else { // when gds > 1
+ } else { // when gds > 1
final double gds = -Math.log((b - p) / k);
- if(Math.log(random.nextDouble()) <= ((k - 1.0) * Math.log(gds))) {
+ if (Math.log(random.nextDouble()) <= ((k - 1.0) * Math.log(gds))) {
return (gds / theta);
}
}
}
- }
- else {
+ } else {
// Step 1. Preparations
final double ss, s, d;
- if(k != -1.0) {
+ if (k != -1.0) {
ss = k - 0.5;
s = Math.sqrt(ss);
d = 5.656854249 - 12.0 * s;
- }
- else {
+ } else {
// For k == -1.0:
ss = 0.0;
s = 0.0;
@@ -465,7 +477,7 @@ public class GammaDistribution implements DistributionWithRandom {
tv2 = 2.0 * random.nextDouble() - 1.0;
tv12 = tv1 * tv1 + tv2 * tv2;
}
- while(tv12 > 1.0);
+ while (tv12 > 1.0);
v1 = tv1;
/* v2 = tv2; */
v12 = tv12;
@@ -480,39 +492,36 @@ public class GammaDistribution implements DistributionWithRandom {
final double t = v1 * Math.sqrt(-2.0 * Math.log(v12) / v12);
final double x = s + 0.5 * t;
final double gds = x * x;
- if(t >= 0.0) {
+ if (t >= 0.0) {
return (gds / theta); // Immediate acceptance
}
// Random uniform
final double un = random.nextDouble();
// Squeeze acceptance
- if(d * un <= t * t * t) {
+ if (d * un <= t * t * t) {
return (gds / theta);
}
- if(k != -1.0) { // Step 4. Set-up for hat case
+ if (k != -1.0) { // Step 4. Set-up for hat case
final double r = 1.0 / k;
q0 = ((((((((q9 * r + q8) * r + q7) * r + q6) * r + q5) * r + q4) * r + q3) * r + q2) * r + q1) * r;
- if(k > 3.686) {
- if(k > 13.022) {
+ if (k > 3.686) {
+ if (k > 13.022) {
b = 1.77;
si = 0.75;
c = 0.1515 / s;
- }
- else {
+ } else {
b = 1.654 + 0.0076 * ss;
si = 1.68 / s + 0.275;
c = 0.062 / s + 0.024;
}
- }
- else {
+ } else {
b = 0.463 + s - 0.178 * ss;
si = 1.235;
c = 0.195 / s - 0.079 + 0.016 * s;
}
- }
- else {
+ } else {
// For k == -1.0:
b = 0.0;
c = 0.0;
@@ -520,24 +529,23 @@ public class GammaDistribution implements DistributionWithRandom {
q0 = 0.0;
}
// Compute v and q
- if(x > 0.0) {
+ if (x > 0.0) {
final double v = t / (s + s);
final double q;
- if(Math.abs(v) > 0.25) {
+ if (Math.abs(v) > 0.25) {
q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
- }
- else {
+ } else {
q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
}
// Quotient acceptance:
- if(Math.log(1.0 - un) <= q) {
+ if (Math.log(1.0 - un) <= q) {
return (gds / theta);
}
}
}
// Double exponential deviate t
- while(true) {
+ while (true) {
double e, u, sign_u, t;
// Retry until t is sufficiently large
do {
@@ -547,30 +555,28 @@ public class GammaDistribution implements DistributionWithRandom {
sign_u = (u > 0) ? 1.0 : -1.0;
t = b + (e * si) * sign_u;
}
- while(t <= -0.71874483771719);
+ while (t <= -0.71874483771719);
// New v(t) and q(t)
final double v = t / (s + s);
final double q;
- if(Math.abs(v) > 0.25) {
+ if (Math.abs(v) > 0.25) {
q = q0 - s * t + 0.25 * t * t + (ss + ss) * Math.log(1.0 + v);
- }
- else {
+ } else {
q = q0 + 0.5 * t * t * ((((((((a9 * v + a8) * v + a7) * v + a6) * v + a5) * v + a4) * v + a3) * v + a2) * v + a1) * v;
}
- if(q <= 0.0) {
+ if (q <= 0.0) {
continue; // retry
}
// Compute w(t)
final double w;
- if(q > 0.5) {
+ if (q > 0.5) {
w = Math.exp(q) - 1.0;
- }
- else {
+ } else {
w = ((((((e7 * q + e6) * q + e5) * q + e4) * q + e3) * q + e2) * q + e1) * q;
}
// Hat acceptance
- if(c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
+ if (c * u * sign_u <= w * Math.exp(e - 0.5 * t * t)) {
final double x = s + 0.5 * t;
return (x * x / theta);
}
@@ -585,7 +591,7 @@ public class GammaDistribution implements DistributionWithRandom {
*
* Reference:
* <p>
- * Algorithm AS 91: The percentage points of the $\chi$ 2 distribution<br />
+ * Algorithm AS 91: The percentage points of the $\chi^2$ distribution<br />
* D.J. Best, D. E. Roberts<br />
* Journal of the Royal Statistical Society. Series C (Applied Statistics)
* </p>
@@ -599,18 +605,18 @@ public class GammaDistribution implements DistributionWithRandom {
protected static double chisquaredProbitApproximation(final double p, double nu, double g) {
final double EPS1 = 1e-2; // Approximation quality
// Sanity checks
- if(Double.isNaN(p) || Double.isNaN(nu)) {
+ if (Double.isNaN(p) || Double.isNaN(nu)) {
return Double.NaN;
}
// Range check
- if(p <= 0) {
+ if (p <= 0) {
return 0;
}
- if(p >= 1) {
+ if (p >= 1) {
return Double.POSITIVE_INFINITY;
}
// Invalid parameters
- if(nu <= 0) {
+ if (nu <= 0) {
return Double.NaN;
}
// Shape of gamma distribution, "XX" in AS 91
@@ -618,7 +624,7 @@ public class GammaDistribution implements DistributionWithRandom {
// For small chi squared values - AS 91
final double logp = Math.log(p);
- if(nu < -1.24 * logp) {
+ if (nu < -1.24 * logp) {
// FIXME: implement and use logGammap1 instead - more stable?
//
// final double lgam1pa = (alpha < 0.5) ? logGammap1(alpha) :
@@ -626,31 +632,29 @@ public class GammaDistribution implements DistributionWithRandom {
// return Math.exp((lgam1pa + logp) / alpha + MathUtil.LOG2);
// This is literal AS 91, above is the GNU R variant.
return Math.pow(p * k * Math.exp(g + k * MathUtil.LOG2), 1 / k);
- }
- else if(nu > 0.32) {
+ } else if (nu > 0.32) {
// Wilson and Hilferty estimate: - AS 91 at 3
final double x = NormalDistribution.quantile(p, 0, 1);
final double p1 = 2. / (9. * nu);
double ch = nu * Math.pow(x * Math.sqrt(p1) + 1 - p1, 3);
// Better approximation for p tending to 1:
- if(ch > 2.2 * nu + 6) {
+ if (ch > 2.2 * nu + 6) {
ch = -2 * (Math.log1p(-p) - (k - 1) * Math.log(0.5 * ch) + g);
}
return ch;
- }
- else {
+ } else {
// nu <= 0.32, AS 91 at 1
final double C7 = 4.67, C8 = 6.66, C9 = 6.73, C10 = 13.32;
final double ag = Math.log1p(-p) + g + (k - 1) * MathUtil.LOG2;
double ch = 0.4;
- while(true) {
+ while (true) {
final double p1 = 1 + ch * (C7 + ch);
final double p2 = ch * (C9 + ch * (C8 + ch));
final double t = -0.5 + (C7 + 2 * ch) / p1 - (C9 + ch * (C10 + 3 * ch)) / p2;
final double delta = (1 - Math.exp(ag + 0.5 * ch) * p2 / p1) / t;
ch -= delta;
- if(Math.abs(delta) > EPS1 * Math.abs(ch)) {
+ if (Math.abs(delta) > EPS1 * Math.abs(ch)) {
return ch;
}
}
@@ -664,7 +668,7 @@ public class GammaDistribution implements DistributionWithRandom {
*
* Reference:
* <p>
- * Algorithm AS 91: The percentage points of the $\chi$^2 distribution<br />
+ * Algorithm AS 91: The percentage points of the $\chi^2$ distribution<br />
* D.J. Best, D. E. Roberts<br />
* Journal of the Royal Statistical Society. Series C (Applied Statistics)
* </p>
@@ -674,34 +678,34 @@ public class GammaDistribution implements DistributionWithRandom {
* @param theta Theta = 1.0/Beta aka. "scaling" parameter
* @return Probit for Gamma distribution
*/
- @Reference(title = "Algorithm AS 91: The percentage points of the $\\chi$^2 distribution", authors = "D.J. Best, D. E. Roberts", booktitle = "Journal of the Royal Statistical Society. Series C (Applied Statistics)")
+ @Reference(title = "Algorithm AS 91: The percentage points of the $\\chi^2$ distribution", authors = "D.J. Best, D. E. Roberts", booktitle = "Journal of the Royal Statistical Society. Series C (Applied Statistics)")
public static double quantile(double p, double k, double theta) {
final double EPS2 = 5e-7; // final precision of AS 91
final int MAXIT = 1000;
// Avoid degenerates
- if(Double.isNaN(p) || Double.isNaN(k) || Double.isNaN(theta)) {
+ if (Double.isNaN(p) || Double.isNaN(k) || Double.isNaN(theta)) {
return Double.NaN;
}
// Range check
- if(p <= 0) {
+ if (p <= 0) {
return 0;
}
- if(p >= 1) {
+ if (p >= 1) {
return Double.POSITIVE_INFINITY;
}
// Shape parameter check
- if(k < 0 || theta <= 0) {
+ if (k < 0 || theta <= 0) {
return Double.NaN;
}
// Corner case - all at 0
- if(k == 0) {
+ if (k == 0) {
return 0.;
}
int max_newton_iterations = 1;
// For small values, ensure some refinement iterations
- if(k < 1e-10) {
+ if (k < 1e-10) {
max_newton_iterations = 7;
}
@@ -713,17 +717,17 @@ public class GammaDistribution implements DistributionWithRandom {
// Second hald of AS 91 follows:
// Refine ChiSquared approximation
chisq: {
- if(Double.isInfinite(ch)) {
+ if (Double.isInfinite(ch)) {
// Cannot refine infinity
max_newton_iterations = 0;
break chisq;
}
- if(ch < EPS2) {
+ if (ch < EPS2) {
// Do not iterate, but refine with newton method
max_newton_iterations = 20;
break chisq;
}
- if(p > 1 - 1e-14 || p < 1e-100) {
+ if (p > 1 - 1e-14 || p < 1e-100) {
// Not in appropriate value range for AS 91
max_newton_iterations = 20;
break chisq;
@@ -732,11 +736,11 @@ public class GammaDistribution implements DistributionWithRandom {
// Phase II: Iteration
final double c = k - 1;
final double ch0 = ch; // backup initial approximation
- for(int i = 1; i <= MAXIT; i++) {
+ for (int i = 1; i <= MAXIT; i++) {
final double q = ch; // previous approximation
final double p1 = 0.5 * ch;
final double p2 = p - regularizedGammaP(k, p1);
- if(Double.isInfinite(p2) || ch <= 0) {
+ if (Double.isInfinite(p2) || ch <= 0) {
ch = ch0;
max_newton_iterations = 27;
break chisq;
@@ -753,11 +757,11 @@ public class GammaDistribution implements DistributionWithRandom {
final double s6 = (120. + c * (346. + 127. * c)) / 5040.;
ch += t * (1 + 0.5 * t * s1 - b * c * (s1 - b * (s2 - b * (s3 - b * (s4 - b * (s5 - b * s6))))));
}
- if(Math.abs(q - ch) < EPS2 * ch) {
+ if (Math.abs(q - ch) < EPS2 * ch) {
break chisq;
}
// Divergence treatment, from GNU R
- if(Math.abs(q - ch) > 0.1 * Math.abs(ch)) {
+ if (Math.abs(q - ch) > 0.1 * Math.abs(ch)) {
ch = ((ch < q) ? 0.9 : 1.1) * q;
}
}
@@ -765,7 +769,7 @@ public class GammaDistribution implements DistributionWithRandom {
// no convergence in MAXIT iterations -- but we add Newton now...
}
double x = 0.5 * ch / theta;
- if(max_newton_iterations > 0) {
+ if (max_newton_iterations > 0) {
// Refine result using final Newton steps.
// TODO: add unit tests that show an improvement! Maybe in logscale only?
x = gammaQuantileNewtonRefinement(Math.log(p), k, theta, max_newton_iterations, x);
@@ -788,33 +792,33 @@ public class GammaDistribution implements DistributionWithRandom {
protected static double gammaQuantileNewtonRefinement(final double logpt, final double k, final double theta, final int maxit, double x) {
final double EPS_N = 1e-15; // Precision threshold
// 0 is not possible, try MIN_NORMAL instead
- if(x <= 0) {
+ if (x <= 0) {
x = Double.MIN_NORMAL;
}
// Current estimation
double logpc = logcdf(x, k, theta);
- if(x == Double.MIN_NORMAL && logpc > logpt * (1. + 1e-7)) {
+ if (x == Double.MIN_NORMAL && logpc > logpt * (1. + 1e-7)) {
return 0.;
}
- if(logpc == Double.NEGATIVE_INFINITY) {
+ if (logpc == Double.NEGATIVE_INFINITY) {
return 0.;
}
// Refine by newton iterations
- for(int i = 0; i < maxit; i++) {
+ for (int i = 0; i < maxit; i++) {
// Error of current approximation
final double logpe = logpc - logpt;
- if(Math.abs(logpe) < Math.abs(EPS_N * logpt)) {
+ if (Math.abs(logpe) < Math.abs(EPS_N * logpt)) {
break;
}
// Step size is controlled by PDF:
final double g = logpdf(x, k, theta);
- if(g == Double.NEGATIVE_INFINITY) {
+ if (g == Double.NEGATIVE_INFINITY) {
break;
}
final double newx = x - logpe * Math.exp(logpc - g);
// New estimate:
logpc = logcdf(newx, k, theta);
- if(Math.abs(logpc - logpt) > Math.abs(logpe) || (i > 0 && Math.abs(logpc - logpt) == Math.abs(logpe))) {
+ if (Math.abs(logpc - logpt) > Math.abs(logpe) || (i > 0 && Math.abs(logpc - logpt) == Math.abs(logpe))) {
// no further improvement
break;
}
@@ -840,21 +844,20 @@ public class GammaDistribution implements DistributionWithRandom {
*/
@Reference(authors = "J. M. Bernando", title = "Algorithm AS 103: Psi (Digamma) Function", booktitle = "Statistical Algorithms")
public static double digamma(double x) {
- if(!(x > 0)) {
+ if (!(x > 0)) {
return Double.NaN;
}
// Method of equation 5:
- if(x <= 1e-5) {
+ if (x <= 1e-5) {
return -EULERS_CONST - 1. / x;
}
// Method of equation 4:
- else if(x > 49.) {
+ else if (x > 49.) {
final double ix2 = 1. / (x * x);
// Partial series expansion
return Math.log(x) - 0.5 / x - ix2 * ((1.0 / 12.) + ix2 * (1.0 / 120. - ix2 / 252.));
// + O(x^8) error
- }
- else {
+ } else {
// Stirling expansion
return digamma(x + 1.) - 1. / x;
}
@@ -869,78 +872,22 @@ public class GammaDistribution implements DistributionWithRandom {
* @return trigamma value
*/
public static double trigamma(double x) {
- if(!(x > 0)) {
+ if (!(x > 0)) {
return Double.NaN;
}
// Method of equation 5:
- if(x <= 1e-5) {
+ if (x <= 1e-5) {
return 1. / (x * x);
}
// Method of equation 4:
- else if(x > 49.) {
+ else if (x > 49.) {
final double ix2 = 1. / (x * x);
// Partial series expansion
return 1 / x - ix2 / 2. + ix2 / x * (1.0 / 6. - ix2 * (1.0 / 30. + ix2 / 42.));
// + O(x^8) error
- }
- else {
+ } else {
// Stirling expansion
return trigamma(x + 1.) - 1. / (x * x);
}
}
-
- /**
- * Mean least squares estimation of Gamma distribution to a set of
- * observations.
- *
- * @param data Data
- * @return Assumed distribution
- */
- public static GammaDistribution estimate(double[] data) {
- return estimate(data, data.length);
- }
-
- /**
- * Mean least squares estimation of Gamma distribution to a set of
- * observations.
- *
- * Reference:
- * <p>
- * Maximum likelihood estimation of the parameters of the gamma distribution
- * and their bias<br />
- * S. C. Choi, R. Wette<br />
- * in: Technometrics
- * </p>
- *
- * @param data Data
- * @param len Length of array
- * @return Assumed distribution
- */
- @Reference(title = "Maximum likelihood estimation of the parameters of the gamma distribution and their bias", authors = "S. C. Choi, R. Wette", booktitle = "Technometrics", url = "http://www.jstor.org/stable/10.2307/1266892")
- public static GammaDistribution estimate(double[] data, int len) {
- double meanx = 0, meanlogx = 0;
- for(int i = 0; i < len; i++) {
- final double logx = Math.log(data[i]);
- final double deltax = data[i] - meanx;
- final double deltalogx = logx - meanlogx;
- meanx += deltax / (i + 1.);
- meanlogx += deltalogx / (i + 1.);
- }
- // Initial approximation
- final double logmeanx = Math.log(meanx);
- final double diff = logmeanx - meanlogx;
- double k = (3 - diff + Math.sqrt((diff - 3) * (diff - 3) + 24 * diff)) / (12 * diff);
-
- // Refine via newton iteration, based on Choi and Wette equation
- while(true) {
- double kdelta = (Math.log(k) - digamma(k) - diff) / (1 / k - trigamma(k));
- if(Math.abs(kdelta) < 1E-8 || Double.isNaN(kdelta)) {
- break;
- }
- k += kdelta;
- }
- // Estimate theta:
- final double theta = k / meanx;
- return new GammaDistribution(k, theta);
- }
-} \ No newline at end of file
+}
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-17 15:02:32 +0000