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 */ 017 package org.apache.commons.math3.distribution; 018 019 import org.apache.commons.math3.exception.NumberIsTooSmallException; 020 import org.apache.commons.math3.exception.util.LocalizedFormats; 021 import org.apache.commons.math3.special.Gamma; 022 import org.apache.commons.math3.special.Beta; 023 import org.apache.commons.math3.util.FastMath; 024 import org.apache.commons.math3.random.RandomGenerator; 025 import org.apache.commons.math3.random.Well19937c; 026 027 /** 028 * Implements the Beta distribution. 029 * 030 * @see <a href="http://en.wikipedia.org/wiki/Beta_distribution">Beta distribution</a> 031 * @version $Id: BetaDistribution.java 1416643 2012-12-03 19:37:14Z tn $ 032 * @since 2.0 (changed to concrete class in 3.0) 033 */ 034 public class BetaDistribution extends AbstractRealDistribution { 035 /** 036 * Default inverse cumulative probability accuracy. 037 * @since 2.1 038 */ 039 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 040 /** Serializable version identifier. */ 041 private static final long serialVersionUID = -1221965979403477668L; 042 /** First shape parameter. */ 043 private final double alpha; 044 /** Second shape parameter. */ 045 private final double beta; 046 /** Normalizing factor used in density computations. 047 * updated whenever alpha or beta are changed. 048 */ 049 private double z; 050 /** Inverse cumulative probability accuracy. */ 051 private final double solverAbsoluteAccuracy; 052 053 /** 054 * Build a new instance. 055 * 056 * @param alpha First shape parameter (must be positive). 057 * @param beta Second shape parameter (must be positive). 058 */ 059 public BetaDistribution(double alpha, double beta) { 060 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 061 } 062 063 /** 064 * Build a new instance. 065 * 066 * @param alpha First shape parameter (must be positive). 067 * @param beta Second shape parameter (must be positive). 068 * @param inverseCumAccuracy Maximum absolute error in inverse 069 * cumulative probability estimates (defaults to 070 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). 071 * @since 2.1 072 */ 073 public BetaDistribution(double alpha, double beta, double inverseCumAccuracy) { 074 this(new Well19937c(), alpha, beta, inverseCumAccuracy); 075 } 076 077 /** 078 * Creates a β distribution. 079 * 080 * @param rng Random number generator. 081 * @param alpha First shape parameter (must be positive). 082 * @param beta Second shape parameter (must be positive). 083 * @param inverseCumAccuracy Maximum absolute error in inverse 084 * cumulative probability estimates (defaults to 085 * {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}). 086 * @since 3.1 087 */ 088 public BetaDistribution(RandomGenerator rng, 089 double alpha, 090 double beta, 091 double inverseCumAccuracy) { 092 super(rng); 093 094 this.alpha = alpha; 095 this.beta = beta; 096 z = Double.NaN; 097 solverAbsoluteAccuracy = inverseCumAccuracy; 098 } 099 100 /** 101 * Access the first shape parameter, {@code alpha}. 102 * 103 * @return the first shape parameter. 104 */ 105 public double getAlpha() { 106 return alpha; 107 } 108 109 /** 110 * Access the second shape parameter, {@code beta}. 111 * 112 * @return the second shape parameter. 113 */ 114 public double getBeta() { 115 return beta; 116 } 117 118 /** Recompute the normalization factor. */ 119 private void recomputeZ() { 120 if (Double.isNaN(z)) { 121 z = Gamma.logGamma(alpha) + Gamma.logGamma(beta) - Gamma.logGamma(alpha + beta); 122 } 123 } 124 125 /** {@inheritDoc} */ 126 public double density(double x) { 127 recomputeZ(); 128 if (x < 0 || x > 1) { 129 return 0; 130 } else if (x == 0) { 131 if (alpha < 1) { 132 throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_0_FOR_SOME_ALPHA, alpha, 1, false); 133 } 134 return 0; 135 } else if (x == 1) { 136 if (beta < 1) { 137 throw new NumberIsTooSmallException(LocalizedFormats.CANNOT_COMPUTE_BETA_DENSITY_AT_1_FOR_SOME_BETA, beta, 1, false); 138 } 139 return 0; 140 } else { 141 double logX = FastMath.log(x); 142 double log1mX = FastMath.log1p(-x); 143 return FastMath.exp((alpha - 1) * logX + (beta - 1) * log1mX - z); 144 } 145 } 146 147 /** {@inheritDoc} */ 148 public double cumulativeProbability(double x) { 149 if (x <= 0) { 150 return 0; 151 } else if (x >= 1) { 152 return 1; 153 } else { 154 return Beta.regularizedBeta(x, alpha, beta); 155 } 156 } 157 158 /** 159 * Return the absolute accuracy setting of the solver used to estimate 160 * inverse cumulative probabilities. 161 * 162 * @return the solver absolute accuracy. 163 * @since 2.1 164 */ 165 @Override 166 protected double getSolverAbsoluteAccuracy() { 167 return solverAbsoluteAccuracy; 168 } 169 170 /** 171 * {@inheritDoc} 172 * 173 * For first shape parameter {@code alpha} and second shape parameter 174 * {@code beta}, the mean is {@code alpha / (alpha + beta)}. 175 */ 176 public double getNumericalMean() { 177 final double a = getAlpha(); 178 return a / (a + getBeta()); 179 } 180 181 /** 182 * {@inheritDoc} 183 * 184 * For first shape parameter {@code alpha} and second shape parameter 185 * {@code beta}, the variance is 186 * {@code (alpha * beta) / [(alpha + beta)^2 * (alpha + beta + 1)]}. 187 */ 188 public double getNumericalVariance() { 189 final double a = getAlpha(); 190 final double b = getBeta(); 191 final double alphabetasum = a + b; 192 return (a * b) / ((alphabetasum * alphabetasum) * (alphabetasum + 1)); 193 } 194 195 /** 196 * {@inheritDoc} 197 * 198 * The lower bound of the support is always 0 no matter the parameters. 199 * 200 * @return lower bound of the support (always 0) 201 */ 202 public double getSupportLowerBound() { 203 return 0; 204 } 205 206 /** 207 * {@inheritDoc} 208 * 209 * The upper bound of the support is always 1 no matter the parameters. 210 * 211 * @return upper bound of the support (always 1) 212 */ 213 public double getSupportUpperBound() { 214 return 1; 215 } 216 217 /** {@inheritDoc} */ 218 public boolean isSupportLowerBoundInclusive() { 219 return false; 220 } 221 222 /** {@inheritDoc} */ 223 public boolean isSupportUpperBoundInclusive() { 224 return false; 225 } 226 227 /** 228 * {@inheritDoc} 229 * 230 * The support of this distribution is connected. 231 * 232 * @return {@code true} 233 */ 234 public boolean isSupportConnected() { 235 return true; 236 } 237 }