BetaDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */
  17. package org.apache.commons.math3.distribution;

  19. import org.apache.commons.math3.exception.NumberIsTooSmallException;
  20. import org.apache.commons.math3.exception.util.LocalizedFormats;
  21. import org.apache.commons.math3.random.RandomGenerator;
  22. import org.apache.commons.math3.random.Well19937c;
  23. import org.apache.commons.math3.special.Beta;
  24. import org.apache.commons.math3.special.Gamma;
  25. import org.apache.commons.math3.util.FastMath;

  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;

  52.     /**
  53.      * Build a new instance.
  54.      * <p>
  55.      * <b>Note:</b> 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.     }

  69.     /**
  70.      * Build a new instance.
  71.      * <p>
  72.      * <b>Note:</b> 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.     }

  90.     /**
  91.      * Creates a &beta; 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.     }

  102.     /**
  103.      * Creates a &beta; 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);

  119.         this.alpha = alpha;
  120.         this.beta = beta;
  121.         z = Double.NaN;
  122.         solverAbsoluteAccuracy = inverseCumAccuracy;
  123.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  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.     }

  249.     /** {@inheritDoc} */
  250.     public boolean isSupportLowerBoundInclusive() {
  251.         return false;
  252.     }

  254.     /** {@inheritDoc} */
  255.     public boolean isSupportUpperBoundInclusive() {
  256.         return false;
  257.     }

  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. }
Created with JaCoCo 0.7.2.201409121644