You are viewing a plain text version of this content. The canonical link for it is here.
Posted to commits@commons.apache.org by er...@apache.org on 2016/03/25 16:19:31 UTC
[46/53] [math] MATH-1158.
MATH-1158.
Method "createSampler" overridden in "GammaDistribution".
Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo
Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/f72b5e65
Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/f72b5e65
Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/f72b5e65
Branch: refs/heads/feature-MATH-1158
Commit: f72b5e65c0c5210cc0be585959ae84fdbc87aced
Parents: adfa016
Author: Gilles <er...@apache.org>
Authored: Sat Mar 12 03:12:31 2016 +0100
Committer: Gilles <er...@apache.org>
Committed: Sat Mar 12 03:12:31 2016 +0100
----------------------------------------------------------------------
.../math4/distribution/GammaDistribution.java | 104 +++++++++++++++++++
1 file changed, 104 insertions(+)
----------------------------------------------------------------------
http://git-wip-us.apache.org/repos/asf/commons-math/blob/f72b5e65/src/main/java/org/apache/commons/math4/distribution/GammaDistribution.java
----------------------------------------------------------------------
diff --git a/src/main/java/org/apache/commons/math4/distribution/GammaDistribution.java b/src/main/java/org/apache/commons/math4/distribution/GammaDistribution.java
index ac231e3..2377df7 100644
--- a/src/main/java/org/apache/commons/math4/distribution/GammaDistribution.java
+++ b/src/main/java/org/apache/commons/math4/distribution/GammaDistribution.java
@@ -20,6 +20,7 @@ import org.apache.commons.math4.exception.NotStrictlyPositiveException;
import org.apache.commons.math4.exception.util.LocalizedFormats;
import org.apache.commons.math4.random.RandomGenerator;
import org.apache.commons.math4.random.Well19937c;
+import org.apache.commons.math4.rng.UniformRandomProvider;
import org.apache.commons.math4.special.Gamma;
import org.apache.commons.math4.util.FastMath;
@@ -152,6 +153,7 @@ public class GammaDistribution extends AbstractRealDistribution {
* {@code scale <= 0}.
* @since 3.3
*/
+ @Deprecated
public GammaDistribution(RandomGenerator rng, double shape, double scale)
throws NotStrictlyPositiveException {
this(rng, shape, scale, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
@@ -170,6 +172,7 @@ public class GammaDistribution extends AbstractRealDistribution {
* {@code scale <= 0}.
* @since 3.1
*/
+ @Deprecated
public GammaDistribution(RandomGenerator rng,
double shape,
double scale,
@@ -419,6 +422,7 @@ public class GammaDistribution extends AbstractRealDistribution {
* @return random value sampled from the Gamma(shape, scale) distribution
*/
@Override
+ @Deprecated
public double sample() {
if (shape < 1) {
// [1]: p. 228, Algorithm GS
@@ -483,4 +487,104 @@ public class GammaDistribution extends AbstractRealDistribution {
}
}
}
+
+ /**
+ * {@inheritDoc}
+ *
+ * <p>
+ * Sampling algorithms:
+ * <ul>
+ * <li>
+ * For {@code 0 < shape < 1}:
+ * <blockquote>
+ * Ahrens, J. H. and Dieter, U.,
+ * <i>Computer methods for sampling from gamma, beta, Poisson and binomial distributions,</i>
+ * Computing, 12, 223-246, 1974.
+ * </blockquote>
+ * </li>
+ * <li>
+ * For {@code shape >= 1}:
+ * <blockquote>
+ * Marsaglia and Tsang, <i>A Simple Method for Generating
+ * Gamma Variables.</i> ACM Transactions on Mathematical Software,
+ * Volume 26 Issue 3, September, 2000.
+ * </blockquote>
+ * </li>
+ * </ul>
+ * </p>
+ */
+ @Override
+ public RealDistribution.Sampler createSampler(final UniformRandomProvider rng) {
+ return new RealDistribution.Sampler() {
+ /** Gaussian sampling. */
+ final RealDistribution.Sampler gaussian = new NormalDistribution().createSampler(rng);
+
+ /** {@inheritDoc} */
+ @Override
+ public double sample() {
+ if (shape < 1) {
+ // [1]: p. 228, Algorithm GS
+
+ while (true) {
+ // Step 1:
+ final double u = rng.nextDouble();
+ final double bGS = 1 + shape / FastMath.E;
+ final double p = bGS * u;
+
+ if (p <= 1) {
+ // Step 2:
+
+ final double x = FastMath.pow(p, 1 / shape);
+ final double u2 = rng.nextDouble();
+
+ if (u2 > FastMath.exp(-x)) {
+ // Reject
+ continue;
+ } else {
+ return scale * x;
+ }
+ } else {
+ // Step 3:
+
+ final double x = -1 * FastMath.log((bGS - p) / shape);
+ final double u2 = rng.nextDouble();
+
+ if (u2 > FastMath.pow(x, shape - 1)) {
+ // Reject
+ continue;
+ } else {
+ return scale * x;
+ }
+ }
+ }
+ }
+
+ // Now shape >= 1
+
+ final double d = shape - 0.333333333333333333;
+ final double c = 1 / (3 * FastMath.sqrt(d));
+
+ while (true) {
+ final double x = gaussian.sample();
+ final double v = (1 + c * x) * (1 + c * x) * (1 + c * x);
+
+ if (v <= 0) {
+ continue;
+ }
+
+ final double x2 = x * x;
+ final double u = rng.nextDouble();
+
+ // Squeeze
+ if (u < 1 - 0.0331 * x2 * x2) {
+ return scale * d * v;
+ }
+
+ if (FastMath.log(u) < 0.5 * x2 + d * (1 - v + FastMath.log(v))) {
+ return scale * d * v;
+ }
+ }
+ }
+ };
+ }
}