I am lost here. Sorry I am still git-newbie. I see no subsequent commit message showing this change has been backed out, but it looks to me like HEAD has it reverted. Is that right? If so, why no subsequent notification? And how can I get and re-apply the changes below so I can do some testing?
Also, re comments else-thread, does the Chi-Square test in RealDistributionAbstractTest succeed regularly and which of the tests below show sensitivity to PRNG seed? Both? One thing to look at is chi-square tests using the TestUtils impl with a lot of bins and healthy sample size. Generally, mistakes in sampling implementations will show systematic problems in certain bins and the TestUtils impl will report where they are. Phil On 12/16/14 1:49 PM, l...@apache.org wrote: > Repository: commons-math > Updated Branches: > refs/heads/MATH-1153 [created] 9ba0a1cad > > > Added Cheng sampling procedure for beta distribution. > > Thanks to Sergei Lebedev for the patch. > > JIRA: MATH-1153 > > Project: http://git-wip-us.apache.org/repos/asf/commons-math/repo > Commit: http://git-wip-us.apache.org/repos/asf/commons-math/commit/9ba0a1ca > Tree: http://git-wip-us.apache.org/repos/asf/commons-math/tree/9ba0a1ca > Diff: http://git-wip-us.apache.org/repos/asf/commons-math/diff/9ba0a1ca > > Branch: refs/heads/MATH-1153 > Commit: 9ba0a1cadf37b0d1c68b01c706a10a042e19e1f6 > Parents: d9b951c > Author: Luc Maisonobe <l...@apache.org> > Authored: Tue Dec 16 21:48:44 2014 +0100 > Committer: Luc Maisonobe <l...@apache.org> > Committed: Tue Dec 16 21:48:44 2014 +0100 > > ---------------------------------------------------------------------- > pom.xml | 3 + > .../math3/distribution/BetaDistribution.java | 139 +++++++++++++++---- > .../distribution/BetaDistributionTest.java | 74 ++++++++++ > 3 files changed, 192 insertions(+), 24 deletions(-) > ---------------------------------------------------------------------- > > > http://git-wip-us.apache.org/repos/asf/commons-math/blob/9ba0a1ca/pom.xml > ---------------------------------------------------------------------- > diff --git a/pom.xml b/pom.xml > index 3ee5db2..667d36e 100644 > --- a/pom.xml > +++ b/pom.xml > @@ -252,6 +252,9 @@ > <name>Piotr Kochanski</name> > </contributor> > <contributor> > + <name>Sergei Lebedev</name> > + </contributor> > + <contributor> > <name>Bob MacCallum</name> > </contributor> > <contributor> > > http://git-wip-us.apache.org/repos/asf/commons-math/blob/9ba0a1ca/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java > ---------------------------------------------------------------------- > diff --git > a/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java > b/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java > index 3f62f64..458fe23 100644 > --- > a/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java > +++ > b/src/main/java/org/apache/commons/math3/distribution/BetaDistribution.java > @@ -23,6 +23,7 @@ import org.apache.commons.math3.random.Well19937c; > import org.apache.commons.math3.special.Beta; > import org.apache.commons.math3.special.Gamma; > import org.apache.commons.math3.util.FastMath; > +import org.apache.commons.math3.util.Precision; > > /** > * Implements the Beta distribution. > @@ -34,10 +35,12 @@ public class BetaDistribution extends > AbstractRealDistribution { > /** > * Default inverse cumulative probability accuracy. > * @since 2.1 > + * @deprecated as of 3.4, this parameter is not used anymore > */ > + @Deprecated > public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; > /** Serializable version identifier. */ > - private static final long serialVersionUID = -1221965979403477668L; > + private static final long serialVersionUID = 20141216L; > /** First shape parameter. */ > private final double alpha; > /** Second shape parameter. */ > @@ -46,8 +49,6 @@ public class BetaDistribution extends > AbstractRealDistribution { > * updated whenever alpha or beta are changed. > */ > private double z; > - /** Inverse cumulative probability accuracy. */ > - private final double solverAbsoluteAccuracy; > > /** > * Build a new instance. > @@ -63,7 +64,7 @@ public class BetaDistribution extends > AbstractRealDistribution { > * @param beta Second shape parameter (must be positive). > */ > public BetaDistribution(double alpha, double beta) { > - this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); > + this(new Well19937c(), alpha, beta); > } > > /** > @@ -82,9 +83,11 @@ public class BetaDistribution extends > AbstractRealDistribution { > * cumulative probability estimates (defaults to > * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). > * @since 2.1 > + * @deprecated as of 3.4, the inverse cumulative accuracy is not used > anymore > */ > + @Deprecated > public BetaDistribution(double alpha, double beta, double > inverseCumAccuracy) { > - this(new Well19937c(), alpha, beta, inverseCumAccuracy); > + this(alpha, beta); > } > > /** > @@ -96,7 +99,11 @@ public class BetaDistribution extends > AbstractRealDistribution { > * @since 3.3 > */ > public BetaDistribution(RandomGenerator rng, double alpha, double beta) { > - this(rng, alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); > + super(rng); > + > + this.alpha = alpha; > + this.beta = beta; > + z = Double.NaN; > } > > /** > @@ -109,17 +116,14 @@ public class BetaDistribution extends > AbstractRealDistribution { > * cumulative probability estimates (defaults to > * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). > * @since 3.1 > + * @deprecated as of 3.4, the inverse cumulative accuracy is not used > anymore > */ > + @Deprecated > public BetaDistribution(RandomGenerator rng, > double alpha, > double beta, > double inverseCumAccuracy) { > - super(rng); > - > - this.alpha = alpha; > - this.beta = beta; > - z = Double.NaN; > - solverAbsoluteAccuracy = inverseCumAccuracy; > + this(rng, alpha, beta); > } > > /** > @@ -188,18 +192,6 @@ public class BetaDistribution extends > AbstractRealDistribution { > } > > /** > - * Return the absolute accuracy setting of the solver used to estimate > - * inverse cumulative probabilities. > - * > - * @return the solver absolute accuracy. > - * @since 2.1 > - */ > - @Override > - protected double getSolverAbsoluteAccuracy() { > - return solverAbsoluteAccuracy; > - } > - > - /** > * {@inheritDoc} > * > * For first shape parameter {@code alpha} and second shape parameter > @@ -266,4 +258,103 @@ public class BetaDistribution extends > AbstractRealDistribution { > public boolean isSupportConnected() { > return true; > } > + > + /** {@inheritDoc} > + * <p> > + * Sampling is performed using Cheng algorithms: > + * </p> > + * <p> > + * R. C. H. Cheng, "Generating beta variates with nonintegral shape > parameters.". > + * Communications of the ACM, 21, 317–322, 1978. > + * </p> > + */ > + @Override > + public double sample() { > + if (FastMath.min(alpha, beta) > 1) { > + return algorithmBB(); > + } else { > + return algorithmBC(); > + } > + } > + > + /** Returns one sample using Cheng's BB algorithm, when both α and > β are greater than 1. > + * @return sampled value > + */ > + private double algorithmBB() { > + final double a = FastMath.min(alpha, beta); > + final double b = FastMath.max(alpha, beta); > + final double newAlpha = a + b; > + final double newBeta = FastMath.sqrt((newAlpha - 2.) / (2. * a * b - > newAlpha)); > + final double gamma = a + 1. / newBeta; > + > + double r; > + double w; > + double t; > + do { > + final double u1 = random.nextDouble(); > + final double u2 = random.nextDouble(); > + final double v = newBeta * FastMath.log(u1 / (1. - u1)); > + w = a * FastMath.exp(v); > + final double newZ = u1 * u1 * u2; > + r = gamma * v - 1.3862944; > + final double s = a + r - w; > + if (s + 2.609438 >= 5 * newZ) { > + break; > + } > + > + t = FastMath.log(newZ); > + if (s >= t) { > + break; > + } > + } while (r + newAlpha * FastMath.log(newAlpha / (b + w)) < t); > + > + w = FastMath.min(w, Double.MAX_VALUE); > + return Precision.equals(a, alpha) ? w / (b + w) : b / (b + w); > + } > + > + /** Returns one sample using Cheng's BC algorithm, when at least one of > α and β is smaller than 1. > + * @return sampled value > + */ > + private double algorithmBC() { > + final double a = FastMath.max(alpha, beta); > + final double b = FastMath.min(alpha, beta); > + final double newAlpha = a + b; > + final double newBeta = 1. / b; > + final double delta = 1. + a - b; > + final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * newBeta > - 0.777778); > + final double k2 = 0.25 + (0.5 + 0.25 / delta) * b; > + > + double w; > + for (;;) { > + final double u1 = random.nextDouble(); > + final double u2 = random.nextDouble(); > + final double y = u1 * u2; > + final double newZ = u1 * y; > + if (u1 < 0.5) { > + if (0.25 * u2 + newZ - y >= k1) { > + continue; > + } > + } else { > + if (newZ <= 0.25) { > + final double v = newBeta * FastMath.log(u1 / (1. - u1)); > + w = a * FastMath.exp(v); > + break; > + } > + > + if (newZ >= k2) { > + continue; > + } > + } > + > + final double v = newBeta * FastMath.log(u1 / (1. - u1)); > + w = a * FastMath.exp(v); > + if (newAlpha * (FastMath.log(newAlpha / (b + w)) + v) - > 1.3862944 >= FastMath.log(newZ)) { > + break; > + } > + } > + > + w = FastMath.min(w, Double.MAX_VALUE); > + return Precision.equals(a, alpha) ? w / (b + w) : b / (b + w); > + } > + > } > > http://git-wip-us.apache.org/repos/asf/commons-math/blob/9ba0a1ca/src/test/java/org/apache/commons/math3/distribution/BetaDistributionTest.java > ---------------------------------------------------------------------- > diff --git > a/src/test/java/org/apache/commons/math3/distribution/BetaDistributionTest.java > > b/src/test/java/org/apache/commons/math3/distribution/BetaDistributionTest.java > index 217ae66..d26c6f0 100644 > --- > a/src/test/java/org/apache/commons/math3/distribution/BetaDistributionTest.java > +++ > b/src/test/java/org/apache/commons/math3/distribution/BetaDistributionTest.java > @@ -16,10 +16,22 @@ > */ > package org.apache.commons.math3.distribution; > > +import java.util.Arrays; > + > +import org.apache.commons.math3.random.RandomGenerator; > +import org.apache.commons.math3.random.Well1024a; > +import org.apache.commons.math3.random.Well19937a; > +import org.apache.commons.math3.stat.StatUtils; > +import org.apache.commons.math3.stat.inference.KolmogorovSmirnovTest; > +import org.apache.commons.math3.stat.inference.TestUtils; > import org.junit.Assert; > import org.junit.Test; > > public class BetaDistributionTest { > + > + static final double[] alphaBetas = {0.1, 1, 10, 100, 1000}; > + static final double epsilon = StatUtils.min(alphaBetas); > + > @Test > public void testCumulative() { > double[] x = new double[]{-0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, > 0.7, 0.8, 0.9, 1.0, 1.1}; > @@ -303,4 +315,66 @@ public class BetaDistributionTest { > Assert.assertEquals(dist.getNumericalMean(), 2.0 / 7.0, tol); > Assert.assertEquals(dist.getNumericalVariance(), 10.0 / (49.0 * > 8.0), tol); > } > + > + @Test > + public void testMomentsSampling() { > + RandomGenerator random = new Well1024a(0x7829862c82fec2dal); > + final int numSamples = 1000; > + for (final double alpha : alphaBetas) { > + for (final double beta : alphaBetas) { > + final BetaDistribution betaDistribution = new > BetaDistribution(random, alpha, beta); > + final double[] observed = new BetaDistribution(alpha, > beta).sample(numSamples); > + Arrays.sort(observed); > + > + final String distribution = String.format("Beta(%.2f, > %.2f)", alpha, beta); > + Assert.assertEquals(String.format("E[%s]", distribution), > + betaDistribution.getNumericalMean(), > + StatUtils.mean(observed), epsilon); > + Assert.assertEquals(String.format("Var[%s]", distribution), > + betaDistribution.getNumericalVariance(), > + StatUtils.variance(observed), epsilon); > + } > + } > + } > + > + @Test > + public void testGoodnessOfFit() { > + RandomGenerator random = new Well19937a(0x237db1db907b089fl); > + final int numSamples = 1000; > + final double level = 0.01; > + for (final double alpha : alphaBetas) { > + for (final double beta : alphaBetas) { > + final BetaDistribution betaDistribution = new > BetaDistribution(random, alpha, beta); > + final double[] observed = > betaDistribution.sample(numSamples); > + Assert.assertFalse("G goodness-of-fit test rejected null at > alpha = " + level, > + gTest(betaDistribution, observed) < > level); > + Assert.assertFalse("KS goodness-of-fit test rejected null at > alpha = " + level, > + new > KolmogorovSmirnovTest(random).kolmogorovSmirnovTest(betaDistribution, > observed) < level); > + } > + } > + } > + > + private double gTest(final RealDistribution expectedDistribution, final > double[] values) { > + final int numBins = values.length / 30; > + final double[] breaks = new double[numBins]; > + for (int b = 0; b < breaks.length; b++) { > + breaks[b] = > expectedDistribution.inverseCumulativeProbability((double) b / numBins); > + } > + > + final long[] observed = new long[numBins]; > + for (final double value : values) { > + int b = 0; > + do { > + b++; > + } while (b < numBins && value >= breaks[b]); > + > + observed[b - 1]++; > + } > + > + final double[] expected = new double[numBins]; > + Arrays.fill(expected, (double) values.length / numBins); > + > + return TestUtils.gTest(expected, observed); > + } > + > } > > --------------------------------------------------------------------- To unsubscribe, e-mail: dev-unsubscr...@commons.apache.org For additional commands, e-mail: dev-h...@commons.apache.org
