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

import org.apache.commons.math3.exception.NumberIsTooSmallException;
import org.apache.commons.math3.exception.util.LocalizedFormats;
import org.apache.commons.math3.random.RandomGenerator;
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.
See Also:
Beta distribution
Since:
2.0 (changed to concrete class in 3.0)/**
 * Implements the Beta distribution.
 *
 * @see <a href="http://en.wikipedia.org/wiki/Beta_distribution">Beta distribution</a>
 * @since 2.0 (changed to concrete class in 3.0)
 */
public class BetaDistribution extends AbstractRealDistribution {
    
Default inverse cumulative probability accuracy.
Since:
2.1/**
     * Default inverse cumulative probability accuracy.
     * @since 2.1
     */
    public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
    
Serializable version identifier. /** Serializable version identifier. */
    private static final long serialVersionUID = -1221965979403477668L;
    
First shape parameter. /** First shape parameter. */
    private final double alpha;
    
Second shape parameter. /** Second shape parameter. */
    private final double beta;
    
Normalizing factor used in density computations.
updated whenever alpha or beta are changed.
/** Normalizing factor used in density computations.
     * updated whenever alpha or beta are changed.
     */
    private double z;
    
Inverse cumulative probability accuracy. /** Inverse cumulative probability accuracy. */
    private final double solverAbsoluteAccuracy;

    
Build a new instance.

Note: this constructor will implicitly create an instance of Well19937c as random generator to be used for sampling only (see sample() and AbstractRealDistribution.sample(int)). In case no sampling is needed for the created distribution, it is advised to pass null as random generator via the appropriate constructors to avoid the additional initialisation overhead. 
Params:
alpha – First shape parameter (must be positive).
beta – Second shape parameter (must be positive)./**
     * Build a new instance.
     * <p>
     * <b>Note:</b> this constructor will implicitly create an instance of
     * {@link Well19937c} as random generator to be used for sampling only (see
     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
     * needed for the created distribution, it is advised to pass {@code null}
     * as random generator via the appropriate constructors to avoid the
     * additional initialisation overhead.
     *
     * @param alpha First shape parameter (must be positive).
     * @param beta Second shape parameter (must be positive).
     */
    public BetaDistribution(double alpha, double beta) {
        this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    }

    
Build a new instance.

Note: this constructor will implicitly create an instance of Well19937c as random generator to be used for sampling only (see sample() and AbstractRealDistribution.sample(int)). In case no sampling is needed for the created distribution, it is advised to pass null as random generator via the appropriate constructors to avoid the additional initialisation overhead. 
Params:
alpha – First shape parameter (must be positive).
beta – Second shape parameter (must be positive).
inverseCumAccuracy – Maximum absolute error in inverse cumulative probability estimates (defaults to DEFAULT_INVERSE_ABSOLUTE_ACCURACY).
Since:
2.1/**
     * Build a new instance.
     * <p>
     * <b>Note:</b> this constructor will implicitly create an instance of
     * {@link Well19937c} as random generator to be used for sampling only (see
     * {@link #sample()} and {@link #sample(int)}). In case no sampling is
     * needed for the created distribution, it is advised to pass {@code null}
     * as random generator via the appropriate constructors to avoid the
     * additional initialisation overhead.
     *
     * @param alpha First shape parameter (must be positive).
     * @param beta Second shape parameter (must be positive).
     * @param inverseCumAccuracy Maximum absolute error in inverse
     * cumulative probability estimates (defaults to
     * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
     * @since 2.1
     */
    public BetaDistribution(double alpha, double beta, double inverseCumAccuracy) {
        this(new Well19937c(), alpha, beta, inverseCumAccuracy);
    }

    
Creates a β distribution.
Params:
rng – Random number generator.
alpha – First shape parameter (must be positive).
beta – Second shape parameter (must be positive).
Since:
3.3/**
     * Creates a &beta; distribution.
     *
     * @param rng Random number generator.
     * @param alpha First shape parameter (must be positive).
     * @param beta Second shape parameter (must be positive).
     * @since 3.3
     */
    public BetaDistribution(RandomGenerator rng, double alpha, double beta) {
        this(rng, alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    }

    
Creates a β distribution.
Params:
rng – Random number generator.
alpha – First shape parameter (must be positive).
beta – Second shape parameter (must be positive).
inverseCumAccuracy – Maximum absolute error in inverse cumulative probability estimates (defaults to DEFAULT_INVERSE_ABSOLUTE_ACCURACY).
Since:
3.1/**
     * Creates a &beta; distribution.
     *
     * @param rng Random number generator.
     * @param alpha First shape parameter (must be positive).
     * @param beta Second shape parameter (must be positive).
     * @param inverseCumAccuracy Maximum absolute error in inverse
     * cumulative probability estimates (defaults to
     * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
     * @since 3.1
     */
    public BetaDistribution(RandomGenerator rng,
                            double alpha,
                            double beta,
                            double inverseCumAccuracy) {
        super(rng);

        this.alpha = alpha;
        this.beta = beta;
        z = Double.NaN;
        solverAbsoluteAccuracy = inverseCumAccuracy;
    }

    
Access the first shape parameter, alpha. 
Returns:
the first shape parameter./**
     * Access the first shape parameter, {@code alpha}.
     *
     * @return the first shape parameter.
     */
    public double getAlpha() {
        return alpha;
    }

    
Access the second shape parameter, beta. 
Returns:
the second shape parameter./**
     * Access the second shape parameter, {@code beta}.
     *
     * @return the second shape parameter.
     */
    public double getBeta() {
        return beta;
    }

    
Recompute the normalization factor. /** Recompute the normalization factor. */
    private void recomputeZ() {
        if (Double.isNaN(z)) {
            z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta);
        }
    }

    
{@inheritDoc} /** {@inheritDoc} */
    public double density(double x) {
        final double logDensity = logDensity(x);
        return logDensity == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logDensity);
    }

    
{@inheritDoc} /** {@inheritDoc} **/
    @Override
    public double logDensity(double x) {
        recomputeZ();
        if (x < 0 || x > 1) {
            return Double.NEGATIVE_INFINITY;
        } else if (x == 0) {
            if (alpha < 1) {
                throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false);
            }
            return Double.NEGATIVE_INFINITY;
        } else if (x == 1) {
            if (beta < 1) {
                throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false);
            }
            return Double.NEGATIVE_INFINITY;
        } else {
            double logX = FastMath.log(x);
            double log1mX = FastMath.log1p(-x);
            return (alpha - 1) * logX + (beta - 1) * log1mX - z;
        }
    }

    
{@inheritDoc} /** {@inheritDoc} */
    public double cumulativeProbability(double x)  {
        if (x <= 0) {
            return 0;
        } else if (x >= 1) {
            return 1;
        } else {
            return Beta.regularizedBeta(x, alpha, beta);
        }
    }

    
Return the absolute accuracy setting of the solver used to estimate
inverse cumulative probabilities.
Returns:
the solver absolute accuracy.
Since:
2.1/**
     * 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 alpha and second shape parameter beta, the mean is alpha / (alpha + beta). /**
     * {@inheritDoc}
     *
     * For first shape parameter {@code alpha} and second shape parameter
     * {@code beta}, the mean is {@code alpha / (alpha + beta)}.
     */
    public double getNumericalMean() {
        final double a = getAlpha();
        return a / (a + getBeta());
    }

    
{@inheritDoc} For first shape parameter alpha and second shape parameter beta, the variance is (alpha * beta) / [(alpha + beta)^2 * (alpha + beta + 1)]. /**
     * {@inheritDoc}
     *
     * For first shape parameter {@code alpha} and second shape parameter
     * {@code beta}, the variance is
     * {@code (alpha * beta) / [(alpha + beta)^2 * (alpha + beta + 1)]}.
     */
    public double getNumericalVariance() {
        final double a = getAlpha();
        final double b = getBeta();
        final double alphabetasum = a + b;
        return (a * b) / ((alphabetasum * alphabetasum) * (alphabetasum + 1));
    }

    
{@inheritDoc}
The lower bound of the support is always 0 no matter the parameters.
Returns:
lower bound of the support (always 0)/**
     * {@inheritDoc}
     *
     * The lower bound of the support is always 0 no matter the parameters.
     *
     * @return lower bound of the support (always 0)
     */
    public double getSupportLowerBound() {
        return 0;
    }

    
{@inheritDoc}
The upper bound of the support is always 1 no matter the parameters.
Returns:
upper bound of the support (always 1)/**
     * {@inheritDoc}
     *
     * The upper bound of the support is always 1 no matter the parameters.
     *
     * @return upper bound of the support (always 1)
     */
    public double getSupportUpperBound() {
        return 1;
    }

    
{@inheritDoc} /** {@inheritDoc} */
    public boolean isSupportLowerBoundInclusive() {
        return false;
    }

    
{@inheritDoc} /** {@inheritDoc} */
    public boolean isSupportUpperBoundInclusive() {
        return false;
    }

    
{@inheritDoc}
The support of this distribution is connected.
Returns:
true/**
     * {@inheritDoc}
     *
     * The support of this distribution is connected.
     *
     * @return {@code true}
     */
    public boolean isSupportConnected() {
        return true;
    }


    
{@inheritDoc}

Sampling is performed using Cheng algorithms:


R. C. H. Cheng, "Generating beta variates with nonintegral shape parameters.".
                Communications of the ACM, 21, 317–322, 1978.

/** {@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() {
        return ChengBetaSampler.sample(random, alpha, beta);
    }

    
Utility class implementing Cheng's algorithms for beta distribution sampling.

R. C. H. Cheng, "Generating beta variates with nonintegral shape parameters.".
                Communications of the ACM, 21, 317–322, 1978.

Since:
3.6/** Utility class implementing Cheng's algorithms for beta distribution sampling.
     * <p>
     * R. C. H. Cheng, "Generating beta variates with nonintegral shape parameters.".
     *                 Communications of the ACM, 21, 317–322, 1978.
     * </p>
     * @since 3.6
     */
    private static final class ChengBetaSampler {

        
Returns one sample using Cheng's sampling algorithm.
Params:
random – random generator to use
alpha – distribution first shape parameter
beta – distribution second shape parameter
Returns:
sampled value/**
         * Returns one sample using Cheng's sampling algorithm.
         * @param random random generator to use
         * @param alpha distribution first shape parameter
         * @param beta distribution second shape parameter
         * @return sampled value
         */
        static double sample(RandomGenerator random, final double alpha, final double beta) {
            final double a = FastMath.min(alpha, beta);
            final double b = FastMath.max(alpha, beta);

            if (a > 1) {
                return algorithmBB(random, alpha, a, b);
            } else {
                return algorithmBC(random, alpha, b, a);
            }
        }

        
Returns one sample using Cheng's BB algorithm, when both α and β are greater than 1.
Params:
random – random generator to use
a0 – distribution first shape parameter (α)
a – min(α, β) where α, β are the two distribution shape parameters
b – max(α, β) where α, β are the two distribution shape parameters
Returns:
sampled value/**
         * Returns one sample using Cheng's BB algorithm, when both &alpha; and &beta; are greater than 1.
         * @param random random generator to use
         * @param a0 distribution first shape parameter (&alpha;)
         * @param a min(&alpha;, &beta;) where &alpha;, &beta; are the two distribution shape parameters
         * @param b max(&alpha;, &beta;) where &alpha;, &beta; are the two distribution shape parameters
         * @return sampled value
         */
        private static double algorithmBB(RandomGenerator random,
                                          final double a0,
                                          final double a,
                                          final double b) {
            final double alpha = a + b;
            final double beta = FastMath.sqrt((alpha - 2.) / (2. * a * b - alpha));
            final double gamma = a + 1. / beta;

            double r;
            double w;
            double t;
            do {
                final double u1 = random.nextDouble();
                final double u2 = random.nextDouble();
                final double v = beta * (FastMath.log(u1) - FastMath.log1p(-u1));
                w = a * FastMath.exp(v);
                final double z = u1 * u1 * u2;
                r = gamma * v - 1.3862944;
                final double s = a + r - w;
                if (s + 2.609438 >= 5 * z) {
                    break;
                }

                t = FastMath.log(z);
                if (s >= t) {
                    break;
                }
            } while (r + alpha * (FastMath.log(alpha) - FastMath.log(b + w)) < t);

            w = FastMath.min(w, Double.MAX_VALUE);
            return Precision.equals(a, a0) ? w / (b + w) : b / (b + w);
        }

        
Returns one sample using Cheng's BC algorithm, when at least one of α and β is smaller than 1.
Params:
random – random generator to use
a0 – distribution first shape parameter (α)
a – max(α, β) where α, β are the two distribution shape parameters
b – min(α, β) where α, β are the two distribution shape parameters
Returns:
sampled value/**
         * Returns one sample using Cheng's BC algorithm, when at least one of &alpha; and &beta; is smaller than 1.
         * @param random random generator to use
         * @param a0 distribution first shape parameter (&alpha;)
         * @param a max(&alpha;, &beta;) where &alpha;, &beta; are the two distribution shape parameters
         * @param b min(&alpha;, &beta;) where &alpha;, &beta; are the two distribution shape parameters
         * @return sampled value
         */
        private static double algorithmBC(RandomGenerator random,
                                          final double a0,
                                          final double a,
                                          final double b) {
            final double alpha = a + b;
            final double beta = 1. / b;
            final double delta = 1. + a - b;
            final double k1 = delta * (0.0138889 + 0.0416667 * b) / (a * beta - 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 z = u1 * y;
                if (u1 < 0.5) {
                    if (0.25 * u2 + z - y >= k1) {
                        continue;
                    }
                } else {
                    if (z <= 0.25) {
                        final double v = beta * (FastMath.log(u1) - FastMath.log1p(-u1));
                        w = a * FastMath.exp(v);
                        break;
                    }

                    if (z >= k2) {
                        continue;
                    }
                }

                final double v = beta * (FastMath.log(u1) - FastMath.log1p(-u1));
                w = a * FastMath.exp(v);
                if (alpha * (FastMath.log(alpha) - FastMath.log(b + w) + v) - 1.3862944 >= FastMath.log(z)) {
                    break;
                }
            }

            w = FastMath.min(w, Double.MAX_VALUE);
            return Precision.equals(a, a0) ? w / (b + w) : b / (b + w);
        }

    }
}