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