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 }