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: 772119 $ $Date: 2009-05-06 05:43:28 -0400 (Wed, 06 May 2009) $ 032 */ 033 public class NormalDistributionImpl extends AbstractContinuousDistribution 034 implements NormalDistribution, Serializable { 035 036 /** Serializable version identifier */ 037 private static final long serialVersionUID = 8589540077390120676L; 038 039 /** &sqrt;(2 π) */ 040 private static final double SQRT2PI = Math.sqrt(2 * Math.PI); 041 042 /** The mean of this distribution. */ 043 private double mean = 0; 044 045 /** The standard deviation of this distribution. */ 046 private double standardDeviation = 1; 047 048 /** 049 * Create a normal distribution using the given mean and standard deviation. 050 * @param mean mean for this distribution 051 * @param sd standard deviation for this distribution 052 */ 053 public NormalDistributionImpl(double mean, double sd){ 054 super(); 055 setMean(mean); 056 setStandardDeviation(sd); 057 } 058 059 /** 060 * Creates normal distribution with the mean equal to zero and standard 061 * deviation equal to one. 062 */ 063 public NormalDistributionImpl(){ 064 this(0.0, 1.0); 065 } 066 067 /** 068 * Access the mean. 069 * @return mean for this distribution 070 */ 071 public double getMean() { 072 return mean; 073 } 074 075 /** 076 * Modify the mean. 077 * @param mean for this distribution 078 */ 079 public void setMean(double mean) { 080 this.mean = mean; 081 } 082 083 /** 084 * Access the standard deviation. 085 * @return standard deviation for this distribution 086 */ 087 public double getStandardDeviation() { 088 return standardDeviation; 089 } 090 091 /** 092 * Modify the standard deviation. 093 * @param sd standard deviation for this distribution 094 * @throws IllegalArgumentException if <code>sd</code> is not positive. 095 */ 096 public void setStandardDeviation(double sd) { 097 if (sd <= 0.0) { 098 throw MathRuntimeException.createIllegalArgumentException( 099 "standard deviation must be positive ({0})", 100 sd); 101 } 102 standardDeviation = sd; 103 } 104 105 /** 106 * Return the probability density for a particular point. 107 * 108 * @param x The point at which the density should be computed. 109 * @return The pdf at point x. 110 */ 111 public double density(Double x) { 112 double x0 = x - getMean(); 113 return Math.exp(-x0 * x0 / (2 * getStandardDeviation() * getStandardDeviation())) / (getStandardDeviation() * SQRT2PI); 114 } 115 116 /** 117 * For this distribution, X, this method returns P(X < <code>x</code>). 118 * @param x the value at which the CDF is evaluated. 119 * @return CDF evaluted at <code>x</code>. 120 * @throws MathException if the algorithm fails to converge; unless 121 * x is more than 20 standard deviations from the mean, in which case the 122 * convergence exception is caught and 0 or 1 is returned. 123 */ 124 public double cumulativeProbability(double x) throws MathException { 125 try { 126 return 0.5 * (1.0 + Erf.erf((x - mean) / 127 (standardDeviation * Math.sqrt(2.0)))); 128 } catch (MaxIterationsExceededException ex) { 129 if (x < (mean - 20 * standardDeviation)) { // JDK 1.5 blows at 38 130 return 0.0d; 131 } else if (x > (mean + 20 * standardDeviation)) { 132 return 1.0d; 133 } else { 134 throw ex; 135 } 136 } 137 } 138 139 /** 140 * For this distribution, X, this method returns the critical point x, such 141 * that P(X < x) = <code>p</code>. 142 * <p> 143 * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and 144 * <code>Double.POSITIVE_INFINITY</code> for p=1.</p> 145 * 146 * @param p the desired probability 147 * @return x, such that P(X < x) = <code>p</code> 148 * @throws MathException if the inverse cumulative probability can not be 149 * computed due to convergence or other numerical errors. 150 * @throws IllegalArgumentException if <code>p</code> is not a valid 151 * probability. 152 */ 153 @Override 154 public double inverseCumulativeProbability(final double p) 155 throws MathException { 156 if (p == 0) { 157 return Double.NEGATIVE_INFINITY; 158 } 159 if (p == 1) { 160 return Double.POSITIVE_INFINITY; 161 } 162 return super.inverseCumulativeProbability(p); 163 } 164 165 /** 166 * Access the domain value lower bound, based on <code>p</code>, used to 167 * bracket a CDF root. This method is used by 168 * {@link #inverseCumulativeProbability(double)} to find critical values. 169 * 170 * @param p the desired probability for the critical value 171 * @return domain value lower bound, i.e. 172 * P(X < <i>lower bound</i>) < <code>p</code> 173 */ 174 @Override 175 protected double getDomainLowerBound(double p) { 176 double ret; 177 178 if (p < .5) { 179 ret = -Double.MAX_VALUE; 180 } else { 181 ret = getMean(); 182 } 183 184 return ret; 185 } 186 187 /** 188 * Access the domain value upper bound, based on <code>p</code>, used to 189 * bracket a CDF root. This method is used by 190 * {@link #inverseCumulativeProbability(double)} to find critical values. 191 * 192 * @param p the desired probability for the critical value 193 * @return domain value upper bound, i.e. 194 * P(X < <i>upper bound</i>) > <code>p</code> 195 */ 196 @Override 197 protected double getDomainUpperBound(double p) { 198 double ret; 199 200 if (p < .5) { 201 ret = getMean(); 202 } else { 203 ret = Double.MAX_VALUE; 204 } 205 206 return ret; 207 } 208 209 /** 210 * Access the initial domain value, based on <code>p</code>, used to 211 * bracket a CDF root. This method is used by 212 * {@link #inverseCumulativeProbability(double)} to find critical values. 213 * 214 * @param p the desired probability for the critical value 215 * @return initial domain value 216 */ 217 @Override 218 protected double getInitialDomain(double p) { 219 double ret; 220 221 if (p < .5) { 222 ret = getMean() - getStandardDeviation(); 223 } else if (p > .5) { 224 ret = getMean() + getStandardDeviation(); 225 } else { 226 ret = getMean(); 227 } 228 229 return ret; 230 } 231 }