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}