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
018package org.apache.commons.math3.distribution;
019
020import org.apache.commons.math3.exception.NotStrictlyPositiveException;
021import org.apache.commons.math3.exception.NumberIsTooLargeException;
022import org.apache.commons.math3.exception.OutOfRangeException;
023import org.apache.commons.math3.exception.util.LocalizedFormats;
024import org.apache.commons.math3.special.Erf;
025import org.apache.commons.math3.util.FastMath;
026import org.apache.commons.math3.random.RandomGenerator;
027import org.apache.commons.math3.random.Well19937c;
028
029/**
030 * Implementation of the normal (gaussian) distribution.
031 *
032 * @see <a href="http://en.wikipedia.org/wiki/Normal_distribution">Normal distribution (Wikipedia)</a>
033 * @see <a href="http://mathworld.wolfram.com/NormalDistribution.html">Normal distribution (MathWorld)</a>
034 * @version $Id: NormalDistribution.java 1535290 2013-10-24 06:58:32Z luc $
035 */
036public class NormalDistribution extends AbstractRealDistribution {
037    /**
038     * Default inverse cumulative probability accuracy.
039     * @since 2.1
040     */
041    public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
042    /** Serializable version identifier. */
043    private static final long serialVersionUID = 8589540077390120676L;
044    /** &radic;(2) */
045    private static final double SQRT2 = FastMath.sqrt(2.0);
046    /** Mean of this distribution. */
047    private final double mean;
048    /** Standard deviation of this distribution. */
049    private final double standardDeviation;
050    /** The value of {@code log(sd) + 0.5*log(2*pi)} stored for faster computation. */
051    private final double logStandardDeviationPlusHalfLog2Pi;
052    /** Inverse cumulative probability accuracy. */
053    private final double solverAbsoluteAccuracy;
054
055    /**
056     * Create a normal distribution with mean equal to zero and standard
057     * deviation equal to one.
058     */
059    public NormalDistribution() {
060        this(0, 1);
061    }
062
063    /**
064     * Create a normal distribution using the given mean and standard deviation.
065     *
066     * @param mean Mean for this distribution.
067     * @param sd Standard deviation for this distribution.
068     * @throws NotStrictlyPositiveException if {@code sd <= 0}.
069     */
070    public NormalDistribution(double mean, double sd)
071        throws NotStrictlyPositiveException {
072        this(mean, sd, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
073    }
074
075    /**
076     * Create a normal distribution using the given mean, standard deviation and
077     * inverse cumulative distribution accuracy.
078     *
079     * @param mean Mean for this distribution.
080     * @param sd Standard deviation for this distribution.
081     * @param inverseCumAccuracy Inverse cumulative probability accuracy.
082     * @throws NotStrictlyPositiveException if {@code sd <= 0}.
083     * @since 2.1
084     */
085    public NormalDistribution(double mean, double sd, double inverseCumAccuracy)
086        throws NotStrictlyPositiveException {
087        this(new Well19937c(), mean, sd, inverseCumAccuracy);
088    }
089
090    /**
091     * Creates a normal distribution.
092     *
093     * @param rng Random number generator.
094     * @param mean Mean for this distribution.
095     * @param sd Standard deviation for this distribution.
096     * @throws NotStrictlyPositiveException if {@code sd <= 0}.
097     * @since 3.3
098     */
099    public NormalDistribution(RandomGenerator rng, double mean, double sd)
100        throws NotStrictlyPositiveException {
101        this(rng, mean, sd, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
102    }
103
104    /**
105     * Creates a normal distribution.
106     *
107     * @param rng Random number generator.
108     * @param mean Mean for this distribution.
109     * @param sd Standard deviation for this distribution.
110     * @param inverseCumAccuracy Inverse cumulative probability accuracy.
111     * @throws NotStrictlyPositiveException if {@code sd <= 0}.
112     * @since 3.1
113     */
114    public NormalDistribution(RandomGenerator rng,
115                              double mean,
116                              double sd,
117                              double inverseCumAccuracy)
118        throws NotStrictlyPositiveException {
119        super(rng);
120
121        if (sd <= 0) {
122            throw new NotStrictlyPositiveException(LocalizedFormats.STANDARD_DEVIATION, sd);
123        }
124
125        this.mean = mean;
126        standardDeviation = sd;
127        logStandardDeviationPlusHalfLog2Pi = FastMath.log(sd) + 0.5 * FastMath.log(2 * FastMath.PI);
128        solverAbsoluteAccuracy = inverseCumAccuracy;
129    }
130
131    /**
132     * Access the mean.
133     *
134     * @return the mean for this distribution.
135     */
136    public double getMean() {
137        return mean;
138    }
139
140    /**
141     * Access the standard deviation.
142     *
143     * @return the standard deviation for this distribution.
144     */
145    public double getStandardDeviation() {
146        return standardDeviation;
147    }
148
149    /** {@inheritDoc} */
150    public double density(double x) {
151        return FastMath.exp(logDensity(x));
152    }
153
154    /** {@inheritDoc} */
155    @Override
156    public double logDensity(double x) {
157        final double x0 = x - mean;
158        final double x1 = x0 / standardDeviation;
159        return -0.5 * x1 * x1 - logStandardDeviationPlusHalfLog2Pi;
160    }
161
162    /**
163     * {@inheritDoc}
164     *
165     * If {@code x} is more than 40 standard deviations from the mean, 0 or 1
166     * is returned, as in these cases the actual value is within
167     * {@code Double.MIN_VALUE} of 0 or 1.
168     */
169    public double cumulativeProbability(double x)  {
170        final double dev = x - mean;
171        if (FastMath.abs(dev) > 40 * standardDeviation) {
172            return dev < 0 ? 0.0d : 1.0d;
173        }
174        return 0.5 * (1 + Erf.erf(dev / (standardDeviation * SQRT2)));
175    }
176
177    /** {@inheritDoc}
178     * @since 3.2
179     */
180    @Override
181    public double inverseCumulativeProbability(final double p) throws OutOfRangeException {
182        if (p < 0.0 || p > 1.0) {
183            throw new OutOfRangeException(p, 0, 1);
184        }
185        return mean + standardDeviation * SQRT2 * Erf.erfInv(2 * p - 1);
186    }
187
188    /**
189     * {@inheritDoc}
190     *
191     * @deprecated See {@link RealDistribution#cumulativeProbability(double,double)}
192     */
193    @Override@Deprecated
194    public double cumulativeProbability(double x0, double x1)
195        throws NumberIsTooLargeException {
196        return probability(x0, x1);
197    }
198
199    /** {@inheritDoc} */
200    @Override
201    public double probability(double x0,
202                              double x1)
203        throws NumberIsTooLargeException {
204        if (x0 > x1) {
205            throw new NumberIsTooLargeException(LocalizedFormats.LOWER_ENDPOINT_ABOVE_UPPER_ENDPOINT,
206                                                x0, x1, true);
207        }
208        final double denom = standardDeviation * SQRT2;
209        final double v0 = (x0 - mean) / denom;
210        final double v1 = (x1 - mean) / denom;
211        return 0.5 * Erf.erf(v0, v1);
212    }
213
214    /** {@inheritDoc} */
215    @Override
216    protected double getSolverAbsoluteAccuracy() {
217        return solverAbsoluteAccuracy;
218    }
219
220    /**
221     * {@inheritDoc}
222     *
223     * For mean parameter {@code mu}, the mean is {@code mu}.
224     */
225    public double getNumericalMean() {
226        return getMean();
227    }
228
229    /**
230     * {@inheritDoc}
231     *
232     * For standard deviation parameter {@code s}, the variance is {@code s^2}.
233     */
234    public double getNumericalVariance() {
235        final double s = getStandardDeviation();
236        return s * s;
237    }
238
239    /**
240     * {@inheritDoc}
241     *
242     * The lower bound of the support is always negative infinity
243     * no matter the parameters.
244     *
245     * @return lower bound of the support (always
246     * {@code Double.NEGATIVE_INFINITY})
247     */
248    public double getSupportLowerBound() {
249        return Double.NEGATIVE_INFINITY;
250    }
251
252    /**
253     * {@inheritDoc}
254     *
255     * The upper bound of the support is always positive infinity
256     * no matter the parameters.
257     *
258     * @return upper bound of the support (always
259     * {@code Double.POSITIVE_INFINITY})
260     */
261    public double getSupportUpperBound() {
262        return Double.POSITIVE_INFINITY;
263    }
264
265    /** {@inheritDoc} */
266    public boolean isSupportLowerBoundInclusive() {
267        return false;
268    }
269
270    /** {@inheritDoc} */
271    public boolean isSupportUpperBoundInclusive() {
272        return false;
273    }
274
275    /**
276     * {@inheritDoc}
277     *
278     * The support of this distribution is connected.
279     *
280     * @return {@code true}
281     */
282    public boolean isSupportConnected() {
283        return true;
284    }
285
286    /** {@inheritDoc} */
287    @Override
288    public double sample()  {
289        return standardDeviation * random.nextGaussian() + mean;
290    }
291}