001/*
002 * Licensed to the Apache Software Foundation (ASF) under one or more
003 * contributor license agreements.  See the NOTICE file distributed with
004 * this work for additional information regarding copyright ownership.
005 * The ASF licenses this file to You under the Apache License, Version 2.0
006 * (the "License"); you may not use this file except in compliance with
007 * the License.  You may obtain a copy of the License at
008 *
009 *      http://www.apache.org/licenses/LICENSE-2.0
010 *
011 * Unless required by applicable law or agreed to in writing, software
012 * distributed under the License is distributed on an "AS IS" BASIS,
013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
014 * See the License for the specific language governing permissions and
015 * limitations under the License.
016 */
017package org.apache.commons.numbers.gamma;
018
019import org.apache.commons.numbers.fraction.ContinuedFraction;
020
021/**
022 * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
023 * Regularized Beta function</a>.
024 * <p>
025 * This class is immutable.
026 * </p>
027 */
028public final class RegularizedBeta {
029    /** Maximum allowed numerical error. */
030    private static final double DEFAULT_EPSILON = 1e-14;
031
032    /** Private constructor. */
033    private RegularizedBeta() {
034        // intentionally empty.
035    }
036
037    /**
038     * Computes the value of the
039     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
040     * regularized beta function</a> I(x, a, b).
041     *
042     * @param x Value.
043     * @param a Parameter {@code a}.
044     * @param b Parameter {@code b}.
045     * @return the regularized beta function I(x, a, b).
046     * @throws ArithmeticException if the algorithm fails to converge.
047     */
048    public static double value(double x,
049                               double a,
050                               double b) {
051        return value(x, a, b, DEFAULT_EPSILON, Integer.MAX_VALUE);
052    }
053
054
055    /**
056     * Computes the value of the
057     * <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
058     * regularized beta function</a> I(x, a, b).
059     *
060     * The implementation of this method is based on:
061     * <ul>
062     *  <li>
063     *   <a href="http://mathworld.wolfram.com/RegularizedBetaFunction.html">
064     *   Regularized Beta Function</a>.
065     *  </li>
066     *  <li>
067     *   <a href="http://functions.wolfram.com/06.21.10.0001.01">
068     *   Regularized Beta Function</a>.
069     *  </li>
070     * </ul>
071     *
072     * @param x the value.
073     * @param a Parameter {@code a}.
074     * @param b Parameter {@code b}.
075     * @param epsilon When the absolute value of the nth item in the
076     * series is less than epsilon the approximation ceases to calculate
077     * further elements in the series.
078     * @param maxIterations Maximum number of "iterations" to complete.
079     * @return the regularized beta function I(x, a, b).
080     * @throws ArithmeticException if the algorithm fails to converge.
081     */
082    public static double value(double x,
083                               final double a,
084                               final double b,
085                               double epsilon,
086                               int maxIterations) {
087        if (Double.isNaN(x) ||
088            Double.isNaN(a) ||
089            Double.isNaN(b) ||
090            x < 0 ||
091            x > 1 ||
092            a <= 0 ||
093            b <= 0) {
094            return Double.NaN;
095        } else if (x > (a + 1) / (2 + b + a) &&
096                   1 - x <= (b + 1) / (2 + b + a)) {
097            return 1 - value(1 - x, b, a, epsilon, maxIterations);
098        } else {
099            final ContinuedFraction fraction = new ContinuedFraction() {
100                /** {@inheritDoc} */
101                @Override
102                protected double getA(int n, double x) {
103                    if (n % 2 == 0) { // even
104                        final double m = n / 2d;
105                        return (m * (b - m) * x) /
106                            ((a + (2 * m) - 1) * (a + (2 * m)));
107                    } else {
108                        final double m = (n - 1d) / 2d;
109                        return -((a + m) * (a + b + m) * x) /
110                            ((a + (2 * m)) * (a + (2 * m) + 1));
111                    }
112                }
113
114                /** {@inheritDoc} */
115                @Override
116                protected double getB(int n, double x) {
117                    return 1;
118                }
119            };
120
121            return Math.exp((a * Math.log(x)) + (b * Math.log1p(-x)) -
122                            Math.log(a) - LogBeta.value(a, b)) /
123                fraction.evaluate(x, epsilon, maxIterations);
124        }
125    }
126}