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.exception.util.LocalizedFormats;
025    import org.apache.commons.math.special.Erf;
026    import org.apache.commons.math.util.FastMath;
027    
028    /**
029     * Default implementation of
030     * {@link org.apache.commons.math.distribution.NormalDistribution}.
031     *
032     * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
033     */
034    public class NormalDistributionImpl extends AbstractContinuousDistribution
035            implements NormalDistribution, Serializable {
036    
037        /**
038         * Default inverse cumulative probability accuracy
039         * @since 2.1
040         */
041        public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
042    
043        /** Serializable version identifier */
044        private static final long serialVersionUID = 8589540077390120676L;
045    
046        /** &sqrt;(2 π) */
047        private static final double SQRT2PI = FastMath.sqrt(2 * FastMath.PI);
048    
049        /** The mean of this distribution. */
050        private double mean = 0;
051    
052        /** The standard deviation of this distribution. */
053        private double standardDeviation = 1;
054    
055        /** Inverse cumulative probability accuracy */
056        private final double solverAbsoluteAccuracy;
057    
058        /**
059         * Create a normal distribution using the given mean and standard deviation.
060         * @param mean mean for this distribution
061         * @param sd standard deviation for this distribution
062         */
063        public NormalDistributionImpl(double mean, double sd){
064            this(mean, sd, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
065        }
066    
067        /**
068         * Create a normal distribution using the given mean, standard deviation and
069         * inverse cumulative distribution accuracy.
070         *
071         * @param mean mean for this distribution
072         * @param sd standard deviation for this distribution
073         * @param inverseCumAccuracy inverse cumulative probability accuracy
074         * @since 2.1
075         */
076        public NormalDistributionImpl(double mean, double sd, double inverseCumAccuracy) {
077            super();
078            setMeanInternal(mean);
079            setStandardDeviationInternal(sd);
080            solverAbsoluteAccuracy = inverseCumAccuracy;
081        }
082    
083        /**
084         * Creates normal distribution with the mean equal to zero and standard
085         * deviation equal to one.
086         */
087        public NormalDistributionImpl(){
088            this(0.0, 1.0);
089        }
090    
091        /**
092         * Access the mean.
093         * @return mean for this distribution
094         */
095        public double getMean() {
096            return mean;
097        }
098    
099        /**
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 &lt; <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 &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         * 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 &lt; <i>lower bound</i>) &lt; <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 &lt; <i>upper bound</i>) &gt; <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    }