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 &lt; <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 &lt; 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 &lt; 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 &lt; <i>lower bound</i>) &lt; <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 &lt; <i>upper bound</i>) &gt; <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    }