1 /* 2 * Licensed to the Apache Software Foundation (ASF) under one or more 3 * contributor license agreements. See the NOTICE file distributed with 4 * this work for additional information regarding copyright ownership. 5 * The ASF licenses this file to You under the Apache License, Version 2.0 6 * (the "License"); you may not use this file except in compliance with 7 * the License. You may obtain a copy of the License at 8 * 9 * http://www.apache.org/licenses/LICENSE-2.0 10 * 11 * Unless required by applicable law or agreed to in writing, software 12 * distributed under the License is distributed on an "AS IS" BASIS, 13 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 14 * See the License for the specific language governing permissions and 15 * limitations under the License. 16 */ 17 package org.apache.commons.math.distribution; 18 19 import java.io.Serializable; 20 21 import org.apache.commons.math.ConvergenceException; 22 import org.apache.commons.math.MathException; 23 import org.apache.commons.math.MathRuntimeException; 24 import org.apache.commons.math.analysis.UnivariateRealFunction; 25 import org.apache.commons.math.analysis.solvers.BrentSolver; 26 import org.apache.commons.math.analysis.solvers.UnivariateRealSolverUtils; 27 import org.apache.commons.math.FunctionEvaluationException; 28 import org.apache.commons.math.exception.util.LocalizedFormats; 29 import org.apache.commons.math.random.RandomDataImpl; 30 import org.apache.commons.math.util.FastMath; 31 32 /** 33 * Base class for continuous distributions. Default implementations are 34 * provided for some of the methods that do not vary from distribution to 35 * distribution. 36 * 37 * @version $Revision: 1073498 $ $Date: 2011-02-22 21:57:26 +0100 (mar. 22 fvr. 2011) $ 38 */ 39 public abstract class AbstractContinuousDistribution 40 extends AbstractDistribution 41 implements ContinuousDistribution, Serializable { 42 43 /** Serializable version identifier */ 44 private static final long serialVersionUID = -38038050983108802L; 45 46 /** 47 * RandomData instance used to generate samples from the distribution 48 * @since 2.2 49 */ 50 protected final RandomDataImpl randomData = new RandomDataImpl(); 51 52 /** 53 * Solver absolute accuracy for inverse cumulative computation 54 * @since 2.1 55 */ 56 private double solverAbsoluteAccuracy = BrentSolver.DEFAULT_ABSOLUTE_ACCURACY; 57 58 /** 59 * Default constructor. 60 */ 61 protected AbstractContinuousDistribution() { 62 super(); 63 } 64 65 /** 66 * Return the probability density for a particular point. 67 * @param x The point at which the density should be computed. 68 * @return The pdf at point x. 69 * @throws MathRuntimeException if the specialized class hasn't implemented this function 70 * @since 2.1 71 */ 72 public double density(double x) throws MathRuntimeException { 73 throw new MathRuntimeException(new UnsupportedOperationException(), 74 LocalizedFormats.NO_DENSITY_FOR_THIS_DISTRIBUTION); 75 } 76 77 /** 78 * For this distribution, X, this method returns the critical point x, such 79 * that P(X < x) = <code>p</code>. 80 * 81 * @param p the desired probability 82 * @return x, such that P(X < x) = <code>p</code> 83 * @throws MathException if the inverse cumulative probability can not be 84 * computed due to convergence or other numerical errors. 85 * @throws IllegalArgumentException if <code>p</code> is not a valid 86 * probability. 87 */ 88 public double inverseCumulativeProbability(final double p) 89 throws MathException { 90 if (p < 0.0 || p > 1.0) { 91 throw MathRuntimeException.createIllegalArgumentException( 92 LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0); 93 } 94 95 // by default, do simple root finding using bracketing and default solver. 96 // subclasses can override if there is a better method. 97 UnivariateRealFunction rootFindingFunction = 98 new UnivariateRealFunction() { 99 public double value(double x) throws FunctionEvaluationException { 100 double ret = Double.NaN; 101 try { 102 ret = cumulativeProbability(x) - p; 103 } catch (MathException ex) { 104 throw new FunctionEvaluationException(x, ex.getSpecificPattern(), ex.getGeneralPattern(), ex.getArguments()); 105 } 106 if (Double.isNaN(ret)) { 107 throw new FunctionEvaluationException(x, LocalizedFormats.CUMULATIVE_PROBABILITY_RETURNED_NAN, x, p); 108 } 109 return ret; 110 } 111 }; 112 113 // Try to bracket root, test domain endpoints if this fails 114 double lowerBound = getDomainLowerBound(p); 115 double upperBound = getDomainUpperBound(p); 116 double[] bracket = null; 117 try { 118 bracket = UnivariateRealSolverUtils.bracket( 119 rootFindingFunction, getInitialDomain(p), 120 lowerBound, upperBound); 121 } catch (ConvergenceException ex) { 122 /* 123 * Check domain endpoints to see if one gives value that is within 124 * the default solver's defaultAbsoluteAccuracy of 0 (will be the 125 * case if density has bounded support and p is 0 or 1). 126 */ 127 if (FastMath.abs(rootFindingFunction.value(lowerBound)) < getSolverAbsoluteAccuracy()) { 128 return lowerBound; 129 } 130 if (FastMath.abs(rootFindingFunction.value(upperBound)) < getSolverAbsoluteAccuracy()) { 131 return upperBound; 132 } 133 // Failed bracket convergence was not because of corner solution 134 throw new MathException(ex); 135 } 136 137 // find root 138 double root = UnivariateRealSolverUtils.solve(rootFindingFunction, 139 // override getSolverAbsoluteAccuracy() to use a Brent solver with 140 // absolute accuracy different from BrentSolver default 141 bracket[0],bracket[1], getSolverAbsoluteAccuracy()); 142 return root; 143 } 144 145 /** 146 * Reseeds the random generator used to generate samples. 147 * 148 * @param seed the new seed 149 * @since 2.2 150 */ 151 public void reseedRandomGenerator(long seed) { 152 randomData.reSeed(seed); 153 } 154 155 /** 156 * Generates a random value sampled from this distribution. The default 157 * implementation uses the 158 * <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method.</a> 159 * 160 * @return random value 161 * @since 2.2 162 * @throws MathException if an error occurs generating the random value 163 */ 164 public double sample() throws MathException { 165 return randomData.nextInversionDeviate(this); 166 } 167 168 /** 169 * Generates a random sample from the distribution. The default implementation 170 * generates the sample by calling {@link #sample()} in a loop. 171 * 172 * @param sampleSize number of random values to generate 173 * @since 2.2 174 * @return an array representing the random sample 175 * @throws MathException if an error occurs generating the sample 176 * @throws IllegalArgumentException if sampleSize is not positive 177 */ 178 public double[] sample(int sampleSize) throws MathException { 179 if (sampleSize <= 0) { 180 MathRuntimeException.createIllegalArgumentException(LocalizedFormats.NOT_POSITIVE_SAMPLE_SIZE, sampleSize); 181 } 182 double[] out = new double[sampleSize]; 183 for (int i = 0; i < sampleSize; i++) { 184 out[i] = sample(); 185 } 186 return out; 187 } 188 189 /** 190 * Access the initial domain value, based on <code>p</code>, used to 191 * bracket a CDF root. This method is used by 192 * {@link #inverseCumulativeProbability(double)} to find critical values. 193 * 194 * @param p the desired probability for the critical value 195 * @return initial domain value 196 */ 197 protected abstract double getInitialDomain(double p); 198 199 /** 200 * Access the domain value lower bound, based on <code>p</code>, used to 201 * bracket a CDF root. This method is used by 202 * {@link #inverseCumulativeProbability(double)} to find critical values. 203 * 204 * @param p the desired probability for the critical value 205 * @return domain value lower bound, i.e. 206 * P(X < <i>lower bound</i>) < <code>p</code> 207 */ 208 protected abstract double getDomainLowerBound(double p); 209 210 /** 211 * Access the domain value upper bound, based on <code>p</code>, used to 212 * bracket a CDF root. This method is used by 213 * {@link #inverseCumulativeProbability(double)} to find critical values. 214 * 215 * @param p the desired probability for the critical value 216 * @return domain value upper bound, i.e. 217 * P(X < <i>upper bound</i>) > <code>p</code> 218 */ 219 protected abstract double getDomainUpperBound(double p); 220 221 /** 222 * Returns the solver absolute accuracy for inverse cumulative computation. 223 * 224 * @return the maximum absolute error in inverse cumulative probability estimates 225 * @since 2.1 226 */ 227 protected double getSolverAbsoluteAccuracy() { 228 return solverAbsoluteAccuracy; 229 } 230 231 } 232