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 18 package org.apache.commons.math.distribution; 19 20 import java.io.Serializable; 21 22 import org.apache.commons.math.MathException; 23 import org.apache.commons.math.MathRuntimeException; 24 import org.apache.commons.math.exception.util.LocalizedFormats; 25 import org.apache.commons.math.special.Erf; 26 import org.apache.commons.math.util.FastMath; 27 28 /** 29 * Default implementation of 30 * {@link org.apache.commons.math.distribution.NormalDistribution}. 31 * 32 * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $ 33 */ 34 public class NormalDistributionImpl extends AbstractContinuousDistribution 35 implements NormalDistribution, Serializable { 36 37 /** 38 * Default inverse cumulative probability accuracy 39 * @since 2.1 40 */ 41 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 42 43 /** Serializable version identifier */ 44 private static final long serialVersionUID = 8589540077390120676L; 45 46 /** &sqrt;(2 π) */ 47 private static final double SQRT2PI = FastMath.sqrt(2 * FastMath.PI); 48 49 /** The mean of this distribution. */ 50 private double mean = 0; 51 52 /** The standard deviation of this distribution. */ 53 private double standardDeviation = 1; 54 55 /** Inverse cumulative probability accuracy */ 56 private final double solverAbsoluteAccuracy; 57 58 /** 59 * Create a normal distribution using the given mean and standard deviation. 60 * @param mean mean for this distribution 61 * @param sd standard deviation for this distribution 62 */ 63 public NormalDistributionImpl(double mean, double sd){ 64 this(mean, sd, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 65 } 66 67 /** 68 * Create a normal distribution using the given mean, standard deviation and 69 * inverse cumulative distribution accuracy. 70 * 71 * @param mean mean for this distribution 72 * @param sd standard deviation for this distribution 73 * @param inverseCumAccuracy inverse cumulative probability accuracy 74 * @since 2.1 75 */ 76 public NormalDistributionImpl(double mean, double sd, double inverseCumAccuracy) { 77 super(); 78 setMeanInternal(mean); 79 setStandardDeviationInternal(sd); 80 solverAbsoluteAccuracy = inverseCumAccuracy; 81 } 82 83 /** 84 * Creates normal distribution with the mean equal to zero and standard 85 * deviation equal to one. 86 */ 87 public NormalDistributionImpl(){ 88 this(0.0, 1.0); 89 } 90 91 /** 92 * Access the mean. 93 * @return mean for this distribution 94 */ 95 public double getMean() { 96 return mean; 97 } 98 99 /** 100 * Modify the mean. 101 * @param mean for this distribution 102 * @deprecated as of 2.1 (class will become immutable in 3.0) 103 */ 104 @Deprecated 105 public void setMean(double mean) { 106 setMeanInternal(mean); 107 } 108 109 /** 110 * Modify the mean. 111 * @param newMean for this distribution 112 */ 113 private void setMeanInternal(double newMean) { 114 this.mean = newMean; 115 } 116 117 /** 118 * Access the standard deviation. 119 * @return standard deviation for this distribution 120 */ 121 public double getStandardDeviation() { 122 return standardDeviation; 123 } 124 125 /** 126 * Modify the standard deviation. 127 * @param sd standard deviation for this distribution 128 * @throws IllegalArgumentException if <code>sd</code> is not positive. 129 * @deprecated as of 2.1 (class will become immutable in 3.0) 130 */ 131 @Deprecated 132 public void setStandardDeviation(double sd) { 133 setStandardDeviationInternal(sd); 134 } 135 136 /** 137 * Modify the standard deviation. 138 * @param sd standard deviation for this distribution 139 * @throws IllegalArgumentException if <code>sd</code> is not positive. 140 */ 141 private void setStandardDeviationInternal(double sd) { 142 if (sd <= 0.0) { 143 throw MathRuntimeException.createIllegalArgumentException( 144 LocalizedFormats.NOT_POSITIVE_STANDARD_DEVIATION, 145 sd); 146 } 147 standardDeviation = sd; 148 } 149 150 /** 151 * Return the probability density for a particular point. 152 * 153 * @param x The point at which the density should be computed. 154 * @return The pdf at point x. 155 * @deprecated 156 */ 157 @Deprecated 158 public double density(Double x) { 159 return density(x.doubleValue()); 160 } 161 162 /** 163 * Returns the probability density for a particular point. 164 * 165 * @param x The point at which the density should be computed. 166 * @return The pdf at point x. 167 * @since 2.1 168 */ 169 @Override 170 public double density(double x) { 171 double x0 = x - mean; 172 return FastMath.exp(-x0 * x0 / (2 * standardDeviation * standardDeviation)) / (standardDeviation * SQRT2PI); 173 } 174 175 /** 176 * For this distribution, X, this method returns P(X < <code>x</code>). 177 * If <code>x</code>is more than 40 standard deviations from the mean, 0 or 1 is returned, 178 * as in these cases the actual value is within <code>Double.MIN_VALUE</code> of 0 or 1. 179 * 180 * @param x the value at which the CDF is evaluated. 181 * @return CDF evaluated at <code>x</code>. 182 * @throws MathException if the algorithm fails to converge 183 */ 184 public double cumulativeProbability(double x) throws MathException { 185 final double dev = x - mean; 186 if (FastMath.abs(dev) > 40 * standardDeviation) { 187 return dev < 0 ? 0.0d : 1.0d; 188 } 189 return 0.5 * (1.0 + Erf.erf(dev / 190 (standardDeviation * FastMath.sqrt(2.0)))); 191 } 192 193 /** 194 * Return the absolute accuracy setting of the solver used to estimate 195 * inverse cumulative probabilities. 196 * 197 * @return the solver absolute accuracy 198 * @since 2.1 199 */ 200 @Override 201 protected double getSolverAbsoluteAccuracy() { 202 return solverAbsoluteAccuracy; 203 } 204 205 /** 206 * For this distribution, X, this method returns the critical point x, such 207 * that P(X < x) = <code>p</code>. 208 * <p> 209 * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and 210 * <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 211 * 212 * @param p the desired probability 213 * @return x, such that P(X < x) = <code>p</code> 214 * @throws MathException if the inverse cumulative probability can not be 215 * computed due to convergence or other numerical errors. 216 * @throws IllegalArgumentException if <code>p</code> is not a valid 217 * probability. 218 */ 219 @Override 220 public double inverseCumulativeProbability(final double p) 221 throws MathException { 222 if (p == 0) { 223 return Double.NEGATIVE_INFINITY; 224 } 225 if (p == 1) { 226 return Double.POSITIVE_INFINITY; 227 } 228 return super.inverseCumulativeProbability(p); 229 } 230 231 /** 232 * Generates a random value sampled from this distribution. 233 * 234 * @return random value 235 * @since 2.2 236 * @throws MathException if an error occurs generating the random value 237 */ 238 @Override 239 public double sample() throws MathException { 240 return randomData.nextGaussian(mean, standardDeviation); 241 } 242 243 /** 244 * Access the domain value lower bound, based on <code>p</code>, used to 245 * bracket a CDF root. This method is used by 246 * {@link #inverseCumulativeProbability(double)} to find critical values. 247 * 248 * @param p the desired probability for the critical value 249 * @return domain value lower bound, i.e. 250 * P(X < <i>lower bound</i>) < <code>p</code> 251 */ 252 @Override 253 protected double getDomainLowerBound(double p) { 254 double ret; 255 256 if (p < .5) { 257 ret = -Double.MAX_VALUE; 258 } else { 259 ret = mean; 260 } 261 262 return ret; 263 } 264 265 /** 266 * Access the domain value upper bound, based on <code>p</code>, used to 267 * bracket a CDF root. This method is used by 268 * {@link #inverseCumulativeProbability(double)} to find critical values. 269 * 270 * @param p the desired probability for the critical value 271 * @return domain value upper bound, i.e. 272 * P(X < <i>upper bound</i>) > <code>p</code> 273 */ 274 @Override 275 protected double getDomainUpperBound(double p) { 276 double ret; 277 278 if (p < .5) { 279 ret = mean; 280 } else { 281 ret = Double.MAX_VALUE; 282 } 283 284 return ret; 285 } 286 287 /** 288 * Access the initial domain value, based on <code>p</code>, used to 289 * bracket a CDF root. This method is used by 290 * {@link #inverseCumulativeProbability(double)} to find critical values. 291 * 292 * @param p the desired probability for the critical value 293 * @return initial domain value 294 */ 295 @Override 296 protected double getInitialDomain(double p) { 297 double ret; 298 299 if (p < .5) { 300 ret = mean - standardDeviation; 301 } else if (p > .5) { 302 ret = mean + standardDeviation; 303 } else { 304 ret = mean; 305 } 306 307 return ret; 308 } 309 310 /** 311 * Returns the lower bound of the support for the distribution. 312 * 313 * The lower bound of the support is always negative infinity 314 * no matter the parameters. 315 * 316 * @return lower bound of the support (always Double.NEGATIVE_INFINITY) 317 * @since 2.2 318 */ 319 public double getSupportLowerBound() { 320 return Double.NEGATIVE_INFINITY; 321 } 322 323 /** 324 * Returns the upper bound of the support for the distribution. 325 * 326 * The upper bound of the support is always positive infinity 327 * no matter the parameters. 328 * 329 * @return upper bound of the support (always Double.POSITIVE_INFINITY) 330 * @since 2.2 331 */ 332 public double getSupportUpperBound() { 333 return Double.POSITIVE_INFINITY; 334 } 335 336 /** 337 * Returns the variance. 338 * 339 * For standard deviation parameter <code>s</code>, 340 * the variance is <code>s^2</code> 341 * 342 * @return the variance 343 * @since 2.2 344 */ 345 public double getNumericalVariance() { 346 final double s = getStandardDeviation(); 347 return s * s; 348 } 349 } 350