    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 &pi;) */
    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    }