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 018 package org.apache.commons.math.distribution; 019 020 import java.io.Serializable; 021 022 import org.apache.commons.math.MathException; 023 import org.apache.commons.math.MathRuntimeException; 024 import org.apache.commons.math.MaxIterationsExceededException; 025 import org.apache.commons.math.special.Erf; 026 027 /** 028 * Default implementation of 029 * {@link org.apache.commons.math.distribution.NormalDistribution}. 030 * 031 * @version $Revision: 925812 $ $Date: 2010-03-21 11:49:31 -0400 (Sun, 21 Mar 2010) $ 032 */ 033 public class NormalDistributionImpl extends AbstractContinuousDistribution 034 implements NormalDistribution, Serializable { 035 036 /** 037 * Default inverse cumulative probability accuracy 038 * @since 2.1 039 */ 040 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9; 041 042 /** Serializable version identifier */ 043 private static final long serialVersionUID = 8589540077390120676L; 044 045 /** &sqrt;(2 π) */ 046 private static final double SQRT2PI = Math.sqrt(2 * Math.PI); 047 048 /** The mean of this distribution. */ 049 private double mean = 0; 050 051 /** The standard deviation of this distribution. */ 052 private double standardDeviation = 1; 053 054 /** Inverse cumulative probability accuracy */ 055 private final double solverAbsoluteAccuracy; 056 057 /** 058 * Create a normal distribution using the given mean and standard deviation. 059 * @param mean mean for this distribution 060 * @param sd standard deviation for this distribution 061 */ 062 public NormalDistributionImpl(double mean, double sd){ 063 this(mean, sd, DEFAULT_INVERSE_ABSOLUTE_ACCURACY); 064 } 065 066 /** 067 * Create a normal distribution using the given mean, standard deviation and 068 * inverse cumulative distribution accuracy. 069 * 070 * @param mean mean for this distribution 071 * @param sd standard deviation for this distribution 072 * @param inverseCumAccuracy inverse cumulative probability accuracy 073 * @since 2.1 074 */ 075 public NormalDistributionImpl(double mean, double sd, double inverseCumAccuracy) { 076 super(); 077 setMeanInternal(mean); 078 setStandardDeviationInternal(sd); 079 solverAbsoluteAccuracy = inverseCumAccuracy; 080 } 081 082 /** 083 * Creates normal distribution with the mean equal to zero and standard 084 * deviation equal to one. 085 */ 086 public NormalDistributionImpl(){ 087 this(0.0, 1.0); 088 } 089 090 /** 091 * Access the mean. 092 * @return mean for this distribution 093 */ 094 public double getMean() { 095 return mean; 096 } 097 098 /** 099 * Modify the mean. 100 * @param mean for this distribution 101 * @deprecated as of 2.1 (class will become immutable in 3.0) 102 */ 103 @Deprecated 104 public void setMean(double mean) { 105 setMeanInternal(mean); 106 } 107 /** 108 * Modify the mean. 109 * @param newMean for this distribution 110 */ 111 private void setMeanInternal(double newMean) { 112 this.mean = newMean; 113 } 114 115 /** 116 * Access the standard deviation. 117 * @return standard deviation for this distribution 118 */ 119 public double getStandardDeviation() { 120 return standardDeviation; 121 } 122 123 /** 124 * Modify the standard deviation. 125 * @param sd standard deviation for this distribution 126 * @throws IllegalArgumentException if <code>sd</code> is not positive. 127 * @deprecated as of 2.1 (class will become immutable in 3.0) 128 */ 129 @Deprecated 130 public void setStandardDeviation(double sd) { 131 setStandardDeviationInternal(sd); 132 } 133 /** 134 * Modify the standard deviation. 135 * @param sd standard deviation for this distribution 136 * @throws IllegalArgumentException if <code>sd</code> is not positive. 137 */ 138 private void setStandardDeviationInternal(double sd) { 139 if (sd <= 0.0) { 140 throw MathRuntimeException.createIllegalArgumentException( 141 "standard deviation must be positive ({0})", 142 sd); 143 } 144 standardDeviation = sd; 145 } 146 147 /** 148 * Return the probability density for a particular point. 149 * 150 * @param x The point at which the density should be computed. 151 * @return The pdf at point x. 152 * @deprecated 153 */ 154 public double density(Double x) { 155 return density(x.doubleValue()); 156 } 157 158 /** 159 * Returns the probability density for a particular point. 160 * 161 * @param x The point at which the density should be computed. 162 * @return The pdf at point x. 163 * @since 2.1 164 */ 165 public double density(double x) { 166 double x0 = x - mean; 167 return Math.exp(-x0 * x0 / (2 * standardDeviation * standardDeviation)) / (standardDeviation * SQRT2PI); 168 } 169 170 /** 171 * For this distribution, X, this method returns P(X < <code>x</code>). 172 * @param x the value at which the CDF is evaluated. 173 * @return CDF evaluted at <code>x</code>. 174 * @throws MathException if the algorithm fails to converge; unless 175 * x is more than 20 standard deviations from the mean, in which case the 176 * convergence exception is caught and 0 or 1 is returned. 177 */ 178 public double cumulativeProbability(double x) throws MathException { 179 try { 180 return 0.5 * (1.0 + Erf.erf((x - mean) / 181 (standardDeviation * Math.sqrt(2.0)))); 182 } catch (MaxIterationsExceededException ex) { 183 if (x < (mean - 20 * standardDeviation)) { // JDK 1.5 blows at 38 184 return 0.0d; 185 } else if (x > (mean + 20 * standardDeviation)) { 186 return 1.0d; 187 } else { 188 throw ex; 189 } 190 } 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 * Access the domain value lower bound, based on <code>p</code>, used to 233 * bracket a CDF root. This method is used by 234 * {@link #inverseCumulativeProbability(double)} to find critical values. 235 * 236 * @param p the desired probability for the critical value 237 * @return domain value lower bound, i.e. 238 * P(X < <i>lower bound</i>) < <code>p</code> 239 */ 240 @Override 241 protected double getDomainLowerBound(double p) { 242 double ret; 243 244 if (p < .5) { 245 ret = -Double.MAX_VALUE; 246 } else { 247 ret = mean; 248 } 249 250 return ret; 251 } 252 253 /** 254 * Access the domain value upper bound, based on <code>p</code>, used to 255 * bracket a CDF root. This method is used by 256 * {@link #inverseCumulativeProbability(double)} to find critical values. 257 * 258 * @param p the desired probability for the critical value 259 * @return domain value upper bound, i.e. 260 * P(X < <i>upper bound</i>) > <code>p</code> 261 */ 262 @Override 263 protected double getDomainUpperBound(double p) { 264 double ret; 265 266 if (p < .5) { 267 ret = mean; 268 } else { 269 ret = Double.MAX_VALUE; 270 } 271 272 return ret; 273 } 274 275 /** 276 * Access the initial domain value, based on <code>p</code>, used to 277 * bracket a CDF root. This method is used by 278 * {@link #inverseCumulativeProbability(double)} to find critical values. 279 * 280 * @param p the desired probability for the critical value 281 * @return initial domain value 282 */ 283 @Override 284 protected double getInitialDomain(double p) { 285 double ret; 286 287 if (p < .5) { 288 ret = mean - standardDeviation; 289 } else if (p > .5) { 290 ret = mean + standardDeviation; 291 } else { 292 ret = mean; 293 } 294 295 return ret; 296 } 297 }