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: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $
032 */
033 public class NormalDistributionImpl extends AbstractContinuousDistribution
034 implements NormalDistribution, Serializable {
035
036 /** Serializable version identifier */
037 private static final long serialVersionUID = 8589540077390120676L;
038
039 /** &sqrt;(2 π) */
040 private static final double SQRT2PI = Math.sqrt(2 * Math.PI);
041
042 /** The mean of this distribution. */
043 private double mean = 0;
044
045 /** The standard deviation of this distribution. */
046 private double standardDeviation = 1;
047
048 /**
049 * Create a normal distribution using the given mean and standard deviation.
050 * @param mean mean for this distribution
051 * @param sd standard deviation for this distribution
052 */
053 public NormalDistributionImpl(double mean, double sd){
054 super();
055 setMean(mean);
056 setStandardDeviation(sd);
057 }
058
059 /**
060 * Creates normal distribution with the mean equal to zero and standard
061 * deviation equal to one.
062 */
063 public NormalDistributionImpl(){
064 this(0.0, 1.0);
065 }
066
067 /**
068 * Access the mean.
069 * @return mean for this distribution
070 */
071 public double getMean() {
072 return mean;
073 }
074
075 /**
076 * Modify the mean.
077 * @param mean for this distribution
078 */
079 public void setMean(double mean) {
080 this.mean = mean;
081 }
082
083 /**
084 * Access the standard deviation.
085 * @return standard deviation for this distribution
086 */
087 public double getStandardDeviation() {
088 return standardDeviation;
089 }
090
091 /**
092 * Modify the standard deviation.
093 * @param sd standard deviation for this distribution
094 * @throws IllegalArgumentException if <code>sd</code> is not positive.
095 */
096 public void setStandardDeviation(double sd) {
097 if (sd <= 0.0) {
098 throw MathRuntimeException.createIllegalArgumentException(
099 "standard deviation must be positive ({0})",
100 sd);
101 }
102 standardDeviation = sd;
103 }
104
105 /**
106 * Return the probability density for a particular point.
107 *
108 * @param x The point at which the density should be computed.
109 * @return The pdf at point x.
110 */
111 public double density(Double x) {
112 double x0 = x - getMean();
113 return Math.exp(-x0 * x0 / (2 * getStandardDeviation() * getStandardDeviation())) / (getStandardDeviation() * SQRT2PI);
114 }
115
116 /**
117 * For this distribution, X, this method returns P(X < <code>x</code>).
118 * @param x the value at which the CDF is evaluated.
119 * @return CDF evaluted at <code>x</code>.
120 * @throws MathException if the algorithm fails to converge; unless
121 * x is more than 20 standard deviations from the mean, in which case the
122 * convergence exception is caught and 0 or 1 is returned.
123 */
124 public double cumulativeProbability(double x) throws MathException {
125 try {
126 return 0.5 * (1.0 + Erf.erf((x - mean) /
127 (standardDeviation * Math.sqrt(2.0))));
128 } catch (MaxIterationsExceededException ex) {
129 if (x < (mean - 20 * standardDeviation)) { // JDK 1.5 blows at 38
130 return 0.0d;
131 } else if (x > (mean + 20 * standardDeviation)) {
132 return 1.0d;
133 } else {
134 throw ex;
135 }
136 }
137 }
138
139 /**
140 * For this distribution, X, this method returns the critical point x, such
141 * that P(X < x) = <code>p</code>.
142 * <p>
143 * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and
144 * <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
145 *
146 * @param p the desired probability
147 * @return x, such that P(X < x) = <code>p</code>
148 * @throws MathException if the inverse cumulative probability can not be
149 * computed due to convergence or other numerical errors.
150 * @throws IllegalArgumentException if <code>p</code> is not a valid
151 * probability.
152 */
153 @Override
154 public double inverseCumulativeProbability(final double p)
155 throws MathException {
156 if (p == 0) {
157 return Double.NEGATIVE_INFINITY;
158 }
159 if (p == 1) {
160 return Double.POSITIVE_INFINITY;
161 }
162 return super.inverseCumulativeProbability(p);
163 }
164
165 /**
166 * Access the domain value lower bound, based on <code>p</code>, used to
167 * bracket a CDF root. This method is used by
168 * {@link #inverseCumulativeProbability(double)} to find critical values.
169 *
170 * @param p the desired probability for the critical value
171 * @return domain value lower bound, i.e.
172 * P(X < <i>lower bound</i>) < <code>p</code>
173 */
174 @Override
175 protected double getDomainLowerBound(double p) {
176 double ret;
177
178 if (p < .5) {
179 ret = -Double.MAX_VALUE;
180 } else {
181 ret = getMean();
182 }
183
184 return ret;
185 }
186
187 /**
188 * Access the domain value upper bound, based on <code>p</code>, used to
189 * bracket a CDF root. This method is used by
190 * {@link #inverseCumulativeProbability(double)} to find critical values.
191 *
192 * @param p the desired probability for the critical value
193 * @return domain value upper bound, i.e.
194 * P(X < <i>upper bound</i>) > <code>p</code>
195 */
196 @Override
197 protected double getDomainUpperBound(double p) {
198 double ret;
199
200 if (p < .5) {
201 ret = getMean();
202 } else {
203 ret = Double.MAX_VALUE;
204 }
205
206 return ret;
207 }
208
209 /**
210 * Access the initial domain value, based on <code>p</code>, used to
211 * bracket a CDF root. This method is used by
212 * {@link #inverseCumulativeProbability(double)} to find critical values.
213 *
214 * @param p the desired probability for the critical value
215 * @return initial domain value
216 */
217 @Override
218 protected double getInitialDomain(double p) {
219 double ret;
220
221 if (p < .5) {
222 ret = getMean() - getStandardDeviation();
223 } else if (p > .5) {
224 ret = getMean() + getStandardDeviation();
225 } else {
226 ret = getMean();
227 }
228
229 return ret;
230 }
231 }