Return-Path: X-Original-To: apmail-commons-notifications-archive@minotaur.apache.org Delivered-To: apmail-commons-notifications-archive@minotaur.apache.org Received: from mail.apache.org (hermes.apache.org [140.211.11.3]) by minotaur.apache.org (Postfix) with SMTP id B698B10420 for ; Sun, 17 May 2015 17:06:03 +0000 (UTC) Received: (qmail 78090 invoked by uid 500); 17 May 2015 17:06:03 -0000 Delivered-To: apmail-commons-notifications-archive@commons.apache.org Received: (qmail 78057 invoked by uid 500); 17 May 2015 17:06:03 -0000 Mailing-List: contact notifications-help@commons.apache.org; run by ezmlm Precedence: bulk List-Help: List-Unsubscribe: List-Post: List-Id: Reply-To: dev@commons.apache.org Delivered-To: mailing list notifications@commons.apache.org Received: (qmail 77967 invoked by uid 99); 17 May 2015 17:06:03 -0000 Received: from eris.apache.org (HELO hades.apache.org) (140.211.11.105) by apache.org (qpsmtpd/0.29) with ESMTP; Sun, 17 May 2015 17:06:03 +0000 Received: from hades.apache.org (localhost [127.0.0.1]) by hades.apache.org (ASF Mail Server at hades.apache.org) with ESMTP id 803B7AC02E3 for ; Sun, 17 May 2015 17:06:03 +0000 (UTC) Content-Type: text/plain; charset="utf-8" MIME-Version: 1.0 Content-Transfer-Encoding: 8bit Subject: svn commit: r951633 [31/49] - in /websites/production/commons/content/proper/commons-math: xref-test/ xref-test/org/apache/commons/math3/ xref-test/org/apache/commons/math3/analysis/ xref-test/org/apache/commons/math3/analysis/differentiation/ xref-tes... Date: Sun, 17 May 2015 17:05:54 -0000 To: notifications@commons.apache.org From: luc@apache.org X-Mailer: svnmailer-1.0.9 Message-Id: <20150517170603.803B7AC02E3@hades.apache.org> Modified: websites/production/commons/content/proper/commons-math/xref/org/apache/commons/math3/distribution/BetaDistribution.html ============================================================================== --- websites/production/commons/content/proper/commons-math/xref/org/apache/commons/math3/distribution/BetaDistribution.html (original) +++ websites/production/commons/content/proper/commons-math/xref/org/apache/commons/math3/distribution/BetaDistribution.html Sun May 17 17:05:50 2015 @@ -31,250 +31,387 @@ 23 import org.apache.commons.math3.special.Beta; 24 import org.apache.commons.math3.special.Gamma; 25 import org.apache.commons.math3.util.FastMath; -26 -27 /** -28 * Implements the Beta distribution. -29 * -30 * @see <a href="http://en.wikipedia.org/wiki/Beta_distribution">Beta distribution</a> -31 * @since 2.0 (changed to concrete class in 3.0) -32 */ -33 public class BetaDistribution extends AbstractRealDistribution { -34 /** -35 * Default inverse cumulative probability accuracy. -36 * @since 2.1 -37 */ -38 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; -39 /** Serializable version identifier. */ -40 private static final long serialVersionUID = -1221965979403477668L; -41 /** First shape parameter. */ -42 private final double alpha; -43 /** Second shape parameter. */ -44 private final double beta; -45 /** Normalizing factor used in density computations. -46 * updated whenever alpha or beta are changed. -47 */ -48 private double z; -49 /** Inverse cumulative probability accuracy. */ -50 private final double solverAbsoluteAccuracy; -51 -52 /** -53 * Build a new instance. -54 * -55 * Note: this constructor will implicitly create an instance of -56 * {@link Well19937c} as random generator to be used for sampling only (see -57 * {@link #sample()} and {@link #sample(int)}). In case no sampling is -58 * needed for the created distribution, it is advised to pass {@code null} -59 * as random generator via the appropriate constructors to avoid the -60 * additional initialisation overhead. -61 * -62 * @param alpha First shape parameter (must be positive). -63 * @param beta Second shape parameter (must be positive). -64 */ -65 public BetaDistribution(double alpha, double beta) { -66 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); -67 } -68 -69 /** -70 * Build a new instance. -71 * -72 * Note: this constructor will implicitly create an instance of -73 * {@link Well19937c} as random generator to be used for sampling only (see -74 * {@link #sample()} and {@link #sample(int)}). In case no sampling is -75 * needed for the created distribution, it is advised to pass {@code null} -76 * as random generator via the appropriate constructors to avoid the -77 * additional initialisation overhead. -78 * -79 * @param alpha First shape parameter (must be positive). -80 * @param beta Second shape parameter (must be positive). -81 * @param inverseCumAccuracy Maximum absolute error in inverse -82 * cumulative probability estimates (defaults to -83 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). -84 * @since 2.1 -85 */ -86 public BetaDistribution(double alpha, double beta, double inverseCumAccuracy) { -87 this(new Well19937c(), alpha, beta, inverseCumAccuracy); -88 } -89 -90 /** -91 * Creates a β distribution. -92 * -93 * @param rng Random number generator. -94 * @param alpha First shape parameter (must be positive). -95 * @param beta Second shape parameter (must be positive). -96 * @since 3.3 -97 */ -98 public BetaDistribution(RandomGenerator rng, double alpha, double beta) { -99 this(rng, alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); -100 } -101 -102 /** -103 * Creates a β distribution. -104 * -105 * @param rng Random number generator. -106 * @param alpha First shape parameter (must be positive). -107 * @param beta Second shape parameter (must be positive). -108 * @param inverseCumAccuracy Maximum absolute error in inverse -109 * cumulative probability estimates (defaults to -110 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). -111 * @since 3.1 -112 */ -113 public BetaDistribution(RandomGenerator rng, -114 double alpha, -115 double beta, -116 double inverseCumAccuracy) { -117 super(rng); -118 -119 this.alpha = alpha; -120 this.beta = beta; -121 z = Double.NaN; -122 solverAbsoluteAccuracy = inverseCumAccuracy; -123 } -124 -125 /** -126 * Access the first shape parameter, {@code alpha}. -127 * -128 * @return the first shape parameter. -129 */ -130 public double getAlpha() { -131 return alpha; -132 } -133 -134 /** -135 * Access the second shape parameter, {@code beta}. -136 * -137 * @return the second shape parameter. -138 */ -139 public double getBeta() { -140 return beta; -141 } -142 -143 /** Recompute the normalization factor. */ -144 private void recomputeZ() { -145 if (Double.isNaN(z)) { -146 z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta); -147 } -148 } -149 -150 /** {@inheritDoc} */ -151 public double density(double x) { -152 final double logDensity = logDensity(x); -153 return logDensity == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logDensity); -154 } -155 -156 /** {@inheritDoc} **/ -157 @Override -158 public double logDensity(double x) { -159 recomputeZ(); -160 if (x < 0 || x > 1) { -161 return Double.NEGATIVE_INFINITY; -162 } else if (x == 0) { -163 if (alpha < 1) { -164 throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false); -165 } -166 return Double.NEGATIVE_INFINITY; -167 } else if (x == 1) { -168 if (beta < 1) { -169 throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false); -170 } -171 return Double.NEGATIVE_INFINITY; -172 } else { -173 double logX = FastMath.log(x); -174 double log1mX = FastMath.log1p(-x); -175 return (alpha - 1) * logX + (beta - 1) * log1mX - z; -176 } -177 } -178 -179 /** {@inheritDoc} */ -180 public double cumulativeProbability(double x) { -181 if (x <= 0) { -182 return 0; -183 } else if (x >= 1) { -184 return 1; -185 } else { -186 return Beta.regularizedBeta(x, alpha, beta); -187 } -188 } -189 -190 /** -191 * Return the absolute accuracy setting of the solver used to estimate -192 * inverse cumulative probabilities. -193 * -194 * @return the solver absolute accuracy. -195 * @since 2.1 -196 */ -197 @Override -198 protected double getSolverAbsoluteAccuracy() { -199 return solverAbsoluteAccuracy; -200 } -201 -202 /** -203 * {@inheritDoc} -204 * -205 * For first shape parameter {@code alpha} and second shape parameter -206 * {@code beta}, the mean is {@code alpha / (alpha + beta)}. -207 */ -208 public double getNumericalMean() { -209 final double a = getAlpha(); -210 return a / (a + getBeta()); -211 } -212 -213 /** -214 * {@inheritDoc} -215 * -216 * For first shape parameter {@code alpha} and second shape parameter -217 * {@code beta}, the variance is -218 * {@code (alpha * beta) / [(alpha + beta)^2 * (alpha + beta + 1)]}. -219 */ -220 public double getNumericalVariance() { -221 final double a = getAlpha(); -222 final double b = getBeta(); -223 final double alphabetasum = a + b; -224 return (a * b) / ((alphabetasum * alphabetasum) * (alphabetasum + 1)); -225 } -226 -227 /** -228 * {@inheritDoc} -229 * -230 * The lower bound of the support is always 0 no matter the parameters. -231 * -232 * @return lower bound of the support (always 0) -233 */ -234 public double getSupportLowerBound() { -235 return 0; -236 } -237 -238 /** -239 * {@inheritDoc} -240 * -241 * The upper bound of the support is always 1 no matter the parameters. -242 * -243 * @return upper bound of the support (always 1) -244 */ -245 public double getSupportUpperBound() { -246 return 1; -247 } -248 -249 /** {@inheritDoc} */ -250 public boolean isSupportLowerBoundInclusive() { -251 return false; -252 } -253 -254 /** {@inheritDoc} */ -255 public boolean isSupportUpperBoundInclusive() { -256 return false; -257 } -258 -259 /** -260 * {@inheritDoc} -261 * -262 * The support of this distribution is connected. -263 * -264 * @return {@code true} -265 */ -266 public boolean isSupportConnected() { -267 return true; -268 } -269 } +26 import org.apache.commons.math3.util.Precision; +27 +28 /** +29 * Implements the Beta distribution. +30 * +31 * @see <a href="http://en.wikipedia.org/wiki/Beta_distribution">Beta distribution</a> +32 * @since 2.0 (changed to concrete class in 3.0) +33 */ +34 public class BetaDistribution extends AbstractRealDistribution { +35 /** +36 * Default inverse cumulative probability accuracy. +37 * @since 2.1 +38 */ +39 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; +40 /** Serializable version identifier. */ +41 private static final long serialVersionUID = -1221965979403477668L; +42 /** First shape parameter. */ +43 private final double alpha; +44 /** Second shape parameter. */ +45 private final double beta; +46 /** Normalizing factor used in density computations. +47 * updated whenever alpha or beta are changed. +48 */ +49 private double z; +50 /** Inverse cumulative probability accuracy. */ +51 private final double solverAbsoluteAccuracy; +52 +53 /** +54 * Build a new instance. +55 * +56 * Note: this constructor will implicitly create an instance of +57 * {@link Well19937c} as random generator to be used for sampling only (see +58 * {@link #sample()} and {@link #sample(int)}). In case no sampling is +59 * needed for the created distribution, it is advised to pass {@code null} +60 * as random generator via the appropriate constructors to avoid the +61 * additional initialisation overhead. +62 * +63 * @param alpha First shape parameter (must be positive). +64 * @param beta Second shape parameter (must be positive). +65 */ +66 public BetaDistribution(double alpha, double beta) { +67 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); +68 } +69 +70 /** +71 * Build a new instance. +72 * +73 * Note: this constructor will implicitly create an instance of +74 * {@link Well19937c} as random generator to be used for sampling only (see +75 * {@link #sample()} and {@link #sample(int)}). In case no sampling is +76 * needed for the created distribution, it is advised to pass {@code null} +77 * as random generator via the appropriate constructors to avoid the +78 * additional initialisation overhead. +79 * +80 * @param alpha First shape parameter (must be positive). +81 * @param beta Second shape parameter (must be positive). +82 * @param inverseCumAccuracy Maximum absolute error in inverse +83 * cumulative probability estimates (defaults to +84 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). +85 * @since 2.1 +86 */ +87 public BetaDistribution(double alpha, double beta, double inverseCumAccuracy) { +88 this(new Well19937c(), alpha, beta, inverseCumAccuracy); +89 } +90 +91 /** +92 * Creates a β distribution. +93 * +94 * @param rng Random number generator. +95 * @param alpha First shape parameter (must be positive). +96 * @param beta Second shape parameter (must be positive). +97 * @since 3.3 +98 */ +99 public BetaDistribution(RandomGenerator rng, double alpha, double beta) { +100 this(rng, alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); +101 } +102 +103 /** +104 * Creates a β distribution. +105 * +106 * @param rng Random number generator. +107 * @param alpha First shape parameter (must be positive). +108 * @param beta Second shape parameter (must be positive). +109 * @param inverseCumAccuracy Maximum absolute error in inverse +110 * cumulative probability estimates (defaults to +111 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). +112 * @since 3.1 +113 */ +114 public BetaDistribution(RandomGenerator rng, +115 double alpha, +116 double beta, +117 double inverseCumAccuracy) { +118 super(rng); +119 +120 this.alpha = alpha; +121 this.beta = beta; +122 z = Double.NaN; +123 solverAbsoluteAccuracy = inverseCumAccuracy; +124 } +125 +126 /** +127 * Access the first shape parameter, {@code alpha}. +128 * +129 * @return the first shape parameter. +130 */ +131 public double getAlpha() { +132 return alpha; +133 } +134 +135 /** +136 * Access the second shape parameter, {@code beta}. +137 * +138 * @return the second shape parameter. +139 */ +140 public double getBeta() { +141 return beta; +142 } +143 +144 /** Recompute the normalization factor. */ +145 private void recomputeZ() { +146 if (Double.isNaN(z)) { +147 z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta); +148 } +149 } +150 +151 /** {@inheritDoc} */ +152 public double density(double x) { +153 final double logDensity = logDensity(x); +154 return logDensity == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logDensity); +155 } +156 +157 /** {@inheritDoc} **/ +158 @Override +159 public double logDensity(double x) { +160 recomputeZ(); +161 if (x < 0 || x > 1) { +162 return Double.NEGATIVE_INFINITY; +163 } else if (x == 0) { +164 if (alpha < 1) { +165 throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false); +166 } +167 return Double.NEGATIVE_INFINITY; +168 } else if (x == 1) { +169 if (beta < 1) { +170 throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false); +171 } +172 return Double.NEGATIVE_INFINITY; +173 } else { +174 double logX = FastMath.log(x); +175 double log1mX = FastMath.log1p(-x); +176 return (alpha - 1) * logX + (beta - 1) * log1mX - z; +177 } +178 } +179 +180 /** {@inheritDoc} */ +181 public double cumulativeProbability(double x) { +182 if (x <= 0) { +183 return 0; +184 } else if (x >= 1) { +185 return 1; +186 } else { +187 return Beta.regularizedBeta(x, alpha, beta); +188 } +189 } +190 +191 /** +192 * Return the absolute accuracy setting of the solver used to estimate +193 * inverse cumulative probabilities. +194 * +195 * @return the solver absolute accuracy. +196 * @since 2.1 +197 */ +198 @Override +199 protected double getSolverAbsoluteAccuracy() { +200 return solverAbsoluteAccuracy; +201 } +202 +203 /** +204 * {@inheritDoc} +205 * +206 * For first shape parameter {@code alpha} and second shape parameter +207 * {@code beta}, the mean is {@code alpha / (alpha + beta)}. +208 */ +209 public double getNumericalMean() { +210 final double a = getAlpha(); +211 return a / (a + getBeta()); +212 } +213 +214 /** +215 * {@inheritDoc} +216 * +217 * For first shape parameter {@code alpha} and second shape parameter +218 * {@code beta}, the variance is +219 * {@code (alpha * beta) / [(alpha + beta)^2 * (alpha + beta + 1)]}. +220 */ +221 public double getNumericalVariance() { +222 final double a = getAlpha(); +223 final double b = getBeta(); +224 final double alphabetasum = a + b; +225 return (a * b) / ((alphabetasum * alphabetasum) * (alphabetasum + 1)); +226 } +227 +228 /** +229 * {@inheritDoc} +230 * +231 * The lower bound of the support is always 0 no matter the parameters. +232 * +233 * @return lower bound of the support (always 0) +234 */ +235 public double getSupportLowerBound() { +236 return 0; +237 } +238 +239 /** +240 * {@inheritDoc} +241 * +242 * The upper bound of the support is always 1 no matter the parameters. +243 * +244 * @return upper bound of the support (always 1) +245 */ +246 public double getSupportUpperBound() { +247 return 1; +248 } +249 +250 /** {@inheritDoc} */ +251 public boolean isSupportLowerBoundInclusive() { +252 return false; +253 } +254 +255 /** {@inheritDoc} */ +256 public boolean isSupportUpperBoundInclusive() { +257 return false; +258 } +259 +260 /** +261 * {@inheritDoc} +262 * +263 * The support of this distribution is connected. +264 * +265 * @return {@code true} +266 */ +267 public boolean isSupportConnected() { +268 return true; +269 } +270 +271 +272 /** {@inheritDoc} +273 * +274 * Sampling is performed using Cheng algorithms: +275 * +276 * +277 * R. C. H. Cheng, "Generating beta variates with nonintegral shape parameters.". +278 * Communications of the ACM, 21, 317Ã¢â¬â322, 1978. +279 * +280 */ +281 @Override +282 public double sample() { +283 return ChengBetaSampler.sample(random, alpha, beta); +284 } +285 +286 /** Utility class implementing Cheng's algorithms for beta distribution sampling. +287 * +288 * R. C. H. Cheng, "Generating beta variates with nonintegral shape parameters.". +289 * Communications of the ACM, 21, 317Ã¢â¬â322, 1978. +290 * +291 * @since 3.6 +292 */ +293 private static final class ChengBetaSampler { +294 +295 /** +296 * Returns one sample using Cheng's sampling algorithm. +297 * @param random random generator to use +298 * @param alpha distribution first shape parameter +299 * @param beta distribution second shape parameter +300 * @return sampled value +301 */ +302 static double sample(RandomGenerator random, final double alpha, final double beta) { +303 final double a = FastMath.min(alpha, beta); +304 final double b = FastMath.max(alpha, beta); +305 +306 if (a > 1) { +307 return algorithmBB(random, alpha, a, b); +308 } else { +309 return algorithmBC(random, alpha, b, a); +310 } +311 } +312 +313 /** +314 * Returns one sample using Cheng's BB algorithm, when both α and β are greater than 1. +315 * @param random random generator to use +316 * @param a0 distribution first shape parameter (α) +317 * @param a min(α, β) where α, β are the two distribution shape parameters +318 * @param b max(α, β) where α, β are the two distribution shape parameters +319 * @return sampled value +320 */ +321 private static double algorithmBB(RandomGenerator random, +322 final double a0, +323 final double a, +324 final double b) { +325 final double alpha = a + b; +326 final double beta = FastMath.sqrt((alpha - 2.) / (2. * a * b - alpha)); +327 final double gamma = a + 1. / beta; +328 +329 double r; +330 double w; +331 double t; +332 do { +333 final double u1 = random.nextDouble(); +334 final double u2 = random.nextDouble(); +335 final double v = beta * (FastMath.log(u1) - FastMath.log1p(-u1)); +336 w = a * FastMath.exp(v); +337 final double z = u1 * u1 * u2; +338 r = gamma * v - 1.3862944; +339 final double s = a + r - w; +340 if (s + 2.609438 >= 5 * z) { +341 break; +342 } +343 +344 t = FastMath.log(z); +345 if (s >= t) { +346 break; +347 } +348 } while (r + alpha * (FastMath.log(alpha) - FastMath.log(b + w)) < t); +349 +350 w = FastMath.min(w, Double.MAX_VALUE); +351 return Precision.equals(a, a0) ? w / (b + w) : b / (b + w); +352 } +353 +354 /** +355 * Returns one sample using Cheng's BC algorithm, when at least one of α and β is smaller than 1. +356 * @param random random generator to use +357 * @param a0 distribution first shape parameter (α) +358 * @param a max(α, β) where α, β are the two distribution shape parameters +359 * @param b min(α, β) where α, β are the two distribution shape parameters +360 * @return sampled value +361 */ +362 private static double algorithmBC(RandomGenerator random, +363 final double a0, +364 final double a, +365 final double b) { +366 final double alpha = a + b; +367 final double beta = 1. / b; +368 final double delta = 1. + a - b; +369 final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * beta - 0.777778); +370 final double k2 = 0.25 + (0.5 + 0.25 / delta) * b; +371 +372 double w; +373 for (;;) { +374 final double u1 = random.nextDouble(); +375 final double u2 = random.nextDouble(); +376 final double y = u1 * u2; +377 final double z = u1 * y; +378 if (u1 < 0.5) { +379 if (0.25 * u2 + z - y >= k1) { +380 continue; +381 } +382 } else { +383 if (z <= 0.25) { +384 final double v = beta * (FastMath.log(u1) - FastMath.log1p(-u1)); +385 w = a * FastMath.exp(v); +386 break; +387 } +388 +389 if (z >= k2) { +390 continue; +391 } +392 } +393 +394 final double v = beta * (FastMath.log(u1) - FastMath.log1p(-u1)); +395 w = a * FastMath.exp(v); +396 if (alpha * (FastMath.log(alpha) - FastMath.log(b + w) + v) - 1.3862944 >= FastMath.log(z)) { +397 break; +398 } +399 } +400 +401 w = FastMath.min(w, Double.MAX_VALUE); +402 return Precision.equals(a, a0) ? w / (b + w) : b / (b + w); +403 } +404 +405 } +406 }