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.statistics.distribution; 018 019import org.apache.commons.numbers.gamma.LogBeta; 020import org.apache.commons.numbers.gamma.RegularizedBeta; 021import org.apache.commons.rng.UniformRandomProvider; 022import org.apache.commons.rng.sampling.distribution.ChengBetaSampler; 023 024/** 025 * Implementation of the beta distribution. 026 * 027 * <p>The probability density function of \( X \) is: 028 * 029 * <p>\[ f(x; \alpha, \beta) = \frac{1}{ B(\alpha, \beta)} x^{\alpha-1} (1-x)^{\beta-1} \] 030 * 031 * <p>for \( \alpha > 0 \), 032 * \( \beta > 0 \), \( x \in [0, 1] \), and 033 * the beta function, \( B \), is a normalization constant: 034 * 035 * <p>\[ B(\alpha, \beta) = \frac{\Gamma(\alpha+\beta)}{\Gamma(\alpha) \Gamma(\beta)} \] 036 * 037 * <p>where \( \Gamma \) is the gamma function. 038 * 039 * <p>\( \alpha \) and \( \beta \) are <em>shape</em> parameters. 040 * 041 * @see <a href="https://en.wikipedia.org/wiki/Beta_distribution">Beta distribution (Wikipedia)</a> 042 * @see <a href="https://mathworld.wolfram.com/BetaDistribution.html">Beta distribution (MathWorld)</a> 043 */ 044public final class BetaDistribution extends AbstractContinuousDistribution { 045 /** First shape parameter. */ 046 private final double alpha; 047 /** Second shape parameter. */ 048 private final double beta; 049 /** Normalizing factor used in log density computations. log(beta(a, b)). */ 050 private final double logBeta; 051 /** Cached value for inverse probability function. */ 052 private final double mean; 053 /** Cached value for inverse probability function. */ 054 private final double variance; 055 056 /** 057 * @param alpha First shape parameter (must be positive). 058 * @param beta Second shape parameter (must be positive). 059 */ 060 private BetaDistribution(double alpha, 061 double beta) { 062 this.alpha = alpha; 063 this.beta = beta; 064 logBeta = LogBeta.value(alpha, beta); 065 final double alphabetasum = alpha + beta; 066 mean = alpha / alphabetasum; 067 variance = (alpha * beta) / ((alphabetasum * alphabetasum) * (alphabetasum + 1)); 068 } 069 070 /** 071 * Creates a beta distribution. 072 * 073 * @param alpha First shape parameter (must be positive). 074 * @param beta Second shape parameter (must be positive). 075 * @return the distribution 076 * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}. 077 */ 078 public static BetaDistribution of(double alpha, 079 double beta) { 080 if (alpha <= 0) { 081 throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, alpha); 082 } 083 if (beta <= 0) { 084 throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE, beta); 085 } 086 return new BetaDistribution(alpha, beta); 087 } 088 089 /** 090 * Gets the first shape parameter of this distribution. 091 * 092 * @return the first shape parameter. 093 */ 094 public double getAlpha() { 095 return alpha; 096 } 097 098 /** 099 * Gets the second shape parameter of this distribution. 100 * 101 * @return the second shape parameter. 102 */ 103 public double getBeta() { 104 return beta; 105 } 106 107 /** {@inheritDoc} 108 * 109 * <p>The density is not defined when {@code x = 0, alpha < 1}, or {@code x = 1, beta < 1}. 110 * In this case the limit of infinity is returned. 111 */ 112 @Override 113 public double density(double x) { 114 if (x < 0 || x > 1) { 115 return 0; 116 } 117 return RegularizedBeta.derivative(x, alpha, beta); 118 } 119 120 /** {@inheritDoc} 121 * 122 * <p>The density is not defined when {@code x = 0, alpha < 1}, or {@code x = 1, beta < 1}. 123 * In this case the limit of infinity is returned. 124 */ 125 @Override 126 public double logDensity(double x) { 127 if (x < 0 || x > 1) { 128 return Double.NEGATIVE_INFINITY; 129 } else if (x == 0) { 130 if (alpha < 1) { 131 // Distribution is not valid when x=0, alpha<1 132 // due to a divide by zero error. 133 // Do not raise an exception and return the limit. 134 return Double.POSITIVE_INFINITY; 135 } 136 // Special case of cancellation: x^(a-1) (1-x)^(b-1) / B(a, b) = 1 / B(a, b) 137 if (alpha == 1) { 138 return -logBeta; 139 } 140 return Double.NEGATIVE_INFINITY; 141 } else if (x == 1) { 142 if (beta < 1) { 143 // Distribution is not valid when x=1, beta<1 144 // due to a divide by zero error. 145 // Do not raise an exception and return the limit. 146 return Double.POSITIVE_INFINITY; 147 } 148 // Special case of cancellation: x^(a-1) (1-x)^(b-1) / B(a, b) = 1 / B(a, b) 149 if (beta == 1) { 150 return -logBeta; 151 } 152 return Double.NEGATIVE_INFINITY; 153 } 154 155 // Log computation 156 final double logX = Math.log(x); 157 final double log1mX = Math.log1p(-x); 158 return (alpha - 1) * logX + (beta - 1) * log1mX - logBeta; 159 } 160 161 /** {@inheritDoc} */ 162 @Override 163 public double cumulativeProbability(double x) { 164 if (x <= 0) { 165 return 0; 166 } else if (x >= 1) { 167 return 1; 168 } else { 169 return RegularizedBeta.value(x, alpha, beta); 170 } 171 } 172 173 /** {@inheritDoc} */ 174 @Override 175 public double survivalProbability(double x) { 176 if (x <= 0) { 177 return 1; 178 } else if (x >= 1) { 179 return 0; 180 } else { 181 return RegularizedBeta.complement(x, alpha, beta); 182 } 183 } 184 185 /** 186 * {@inheritDoc} 187 * 188 * <p>For first shape parameter \( \alpha \) and second shape parameter 189 * \( \beta \), the mean is: 190 * 191 * <p>\[ \frac{\alpha}{\alpha + \beta} \] 192 */ 193 @Override 194 public double getMean() { 195 return mean; 196 } 197 198 /** 199 * {@inheritDoc} 200 * 201 * <p>For first shape parameter \( \alpha \) and second shape parameter 202 * \( \beta \), the variance is: 203 * 204 * <p>\[ \frac{\alpha \beta}{(\alpha + \beta)^2 (\alpha + \beta + 1)} \] 205 */ 206 @Override 207 public double getVariance() { 208 return variance; 209 } 210 211 /** 212 * {@inheritDoc} 213 * 214 * <p>The lower bound of the support is always 0. 215 * 216 * @return 0. 217 */ 218 @Override 219 public double getSupportLowerBound() { 220 return 0; 221 } 222 223 /** 224 * {@inheritDoc} 225 * 226 * <p>The upper bound of the support is always 1. 227 * 228 * @return 1. 229 */ 230 @Override 231 public double getSupportUpperBound() { 232 return 1; 233 } 234 235 /** {@inheritDoc} */ 236 @Override 237 public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) { 238 // Beta distribution sampler. 239 return ChengBetaSampler.of(rng, alpha, beta)::sample; 240 } 241}