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.OutOfRangeException;
021 import org.apache.commons.math3.exception.NotStrictlyPositiveException;
022 import org.apache.commons.math3.exception.util.LocalizedFormats;
023 import org.apache.commons.math3.special.Gamma;
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 Weibull distribution. This implementation uses the
030 * two parameter form of the distribution defined by
031 * <a href="http://mathworld.wolfram.com/WeibullDistribution.html">
032 * Weibull Distribution</a>, equations (1) and (2).
033 *
034 * @see <a href="http://en.wikipedia.org/wiki/Weibull_distribution">Weibull distribution (Wikipedia)</a>
035 * @see <a href="http://mathworld.wolfram.com/WeibullDistribution.html">Weibull distribution (MathWorld)</a>
036 * @since 1.1 (changed to concrete class in 3.0)
037 * @version $Id: WeibullDistribution.java 1416643 2012-12-03 19:37:14Z tn $
038 */
039 public class WeibullDistribution extends AbstractRealDistribution {
040 /**
041 * Default inverse cumulative probability accuracy.
042 * @since 2.1
043 */
044 public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
045 /** Serializable version identifier. */
046 private static final long serialVersionUID = 8589540077390120676L;
047 /** The shape parameter. */
048 private final double shape;
049 /** The scale parameter. */
050 private final double scale;
051 /** Inverse cumulative probability accuracy. */
052 private final double solverAbsoluteAccuracy;
053 /** Cached numerical mean */
054 private double numericalMean = Double.NaN;
055 /** Whether or not the numerical mean has been calculated */
056 private boolean numericalMeanIsCalculated = false;
057 /** Cached numerical variance */
058 private double numericalVariance = Double.NaN;
059 /** Whether or not the numerical variance has been calculated */
060 private boolean numericalVarianceIsCalculated = false;
061
062 /**
063 * Create a Weibull distribution with the given shape and scale and a
064 * location equal to zero.
065 *
066 * @param alpha Shape parameter.
067 * @param beta Scale parameter.
068 * @throws NotStrictlyPositiveException if {@code alpha <= 0} or
069 * {@code beta <= 0}.
070 */
071 public WeibullDistribution(double alpha, double beta)
072 throws NotStrictlyPositiveException {
073 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
074 }
075
076 /**
077 * Create a Weibull distribution with the given shape, scale and inverse
078 * cumulative probability accuracy and a location equal to zero.
079 *
080 * @param alpha Shape parameter.
081 * @param beta Scale parameter.
082 * @param inverseCumAccuracy Maximum absolute error in inverse
083 * cumulative probability estimates
084 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
085 * @throws NotStrictlyPositiveException if {@code alpha <= 0} or
086 * {@code beta <= 0}.
087 * @since 2.1
088 */
089 public WeibullDistribution(double alpha, double beta,
090 double inverseCumAccuracy) {
091 this(new Well19937c(), alpha, beta, inverseCumAccuracy);
092 }
093
094 /**
095 * Creates a Weibull distribution.
096 *
097 * @param rng Random number generator.
098 * @param alpha Shape parameter.
099 * @param beta Scale parameter.
100 * @param inverseCumAccuracy Maximum absolute error in inverse
101 * cumulative probability estimates
102 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
103 * @throws NotStrictlyPositiveException if {@code alpha <= 0} or
104 * {@code beta <= 0}.
105 * @since 3.1
106 */
107 public WeibullDistribution(RandomGenerator rng,
108 double alpha,
109 double beta,
110 double inverseCumAccuracy)
111 throws NotStrictlyPositiveException {
112 super(rng);
113
114 if (alpha <= 0) {
115 throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE,
116 alpha);
117 }
118 if (beta <= 0) {
119 throw new NotStrictlyPositiveException(LocalizedFormats.SCALE,
120 beta);
121 }
122 scale = beta;
123 shape = alpha;
124 solverAbsoluteAccuracy = inverseCumAccuracy;
125 }
126
127 /**
128 * Access the shape parameter, {@code alpha}.
129 *
130 * @return the shape parameter, {@code alpha}.
131 */
132 public double getShape() {
133 return shape;
134 }
135
136 /**
137 * Access the scale parameter, {@code beta}.
138 *
139 * @return the scale parameter, {@code beta}.
140 */
141 public double getScale() {
142 return scale;
143 }
144
145 /** {@inheritDoc} */
146 public double density(double x) {
147 if (x < 0) {
148 return 0;
149 }
150
151 final double xscale = x / scale;
152 final double xscalepow = FastMath.pow(xscale, shape - 1);
153
154 /*
155 * FastMath.pow(x / scale, shape) =
156 * FastMath.pow(xscale, shape) =
157 * FastMath.pow(xscale, shape - 1) * xscale
158 */
159 final double xscalepowshape = xscalepow * xscale;
160
161 return (shape / scale) * xscalepow * FastMath.exp(-xscalepowshape);
162 }
163
164 /** {@inheritDoc} */
165 public double cumulativeProbability(double x) {
166 double ret;
167 if (x <= 0.0) {
168 ret = 0.0;
169 } else {
170 ret = 1.0 - FastMath.exp(-FastMath.pow(x / scale, shape));
171 }
172 return ret;
173 }
174
175 /**
176 * {@inheritDoc}
177 *
178 * Returns {@code 0} when {@code p == 0} and
179 * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
180 */
181 @Override
182 public double inverseCumulativeProbability(double p) {
183 double ret;
184 if (p < 0.0 || p > 1.0) {
185 throw new OutOfRangeException(p, 0.0, 1.0);
186 } else if (p == 0) {
187 ret = 0.0;
188 } else if (p == 1) {
189 ret = Double.POSITIVE_INFINITY;
190 } else {
191 ret = scale * FastMath.pow(-FastMath.log(1.0 - p), 1.0 / shape);
192 }
193 return ret;
194 }
195
196 /**
197 * Return the absolute accuracy setting of the solver used to estimate
198 * inverse cumulative probabilities.
199 *
200 * @return the solver absolute accuracy.
201 * @since 2.1
202 */
203 @Override
204 protected double getSolverAbsoluteAccuracy() {
205 return solverAbsoluteAccuracy;
206 }
207
208 /**
209 * {@inheritDoc}
210 *
211 * The mean is {@code scale * Gamma(1 + (1 / shape))}, where {@code Gamma()}
212 * is the Gamma-function.
213 */
214 public double getNumericalMean() {
215 if (!numericalMeanIsCalculated) {
216 numericalMean = calculateNumericalMean();
217 numericalMeanIsCalculated = true;
218 }
219 return numericalMean;
220 }
221
222 /**
223 * used by {@link #getNumericalMean()}
224 *
225 * @return the mean of this distribution
226 */
227 protected double calculateNumericalMean() {
228 final double sh = getShape();
229 final double sc = getScale();
230
231 return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
232 }
233
234 /**
235 * {@inheritDoc}
236 *
237 * The variance is {@code scale^2 * Gamma(1 + (2 / shape)) - mean^2}
238 * where {@code Gamma()} is the Gamma-function.
239 */
240 public double getNumericalVariance() {
241 if (!numericalVarianceIsCalculated) {
242 numericalVariance = calculateNumericalVariance();
243 numericalVarianceIsCalculated = true;
244 }
245 return numericalVariance;
246 }
247
248 /**
249 * used by {@link #getNumericalVariance()}
250 *
251 * @return the variance of this distribution
252 */
253 protected double calculateNumericalVariance() {
254 final double sh = getShape();
255 final double sc = getScale();
256 final double mn = getNumericalMean();
257
258 return (sc * sc) * FastMath.exp(Gamma.logGamma(1 + (2 / sh))) -
259 (mn * mn);
260 }
261
262 /**
263 * {@inheritDoc}
264 *
265 * The lower bound of the support is always 0 no matter the parameters.
266 *
267 * @return lower bound of the support (always 0)
268 */
269 public double getSupportLowerBound() {
270 return 0;
271 }
272
273 /**
274 * {@inheritDoc}
275 *
276 * The upper bound of the support is always positive infinity
277 * no matter the parameters.
278 *
279 * @return upper bound of the support (always
280 * {@code Double.POSITIVE_INFINITY})
281 */
282 public double getSupportUpperBound() {
283 return Double.POSITIVE_INFINITY;
284 }
285
286 /** {@inheritDoc} */
287 public boolean isSupportLowerBoundInclusive() {
288 return true;
289 }
290
291 /** {@inheritDoc} */
292 public boolean isSupportUpperBoundInclusive() {
293 return false;
294 }
295
296 /**
297 * {@inheritDoc}
298 *
299 * The support of this distribution is connected.
300 *
301 * @return {@code true}
302 */
303 public boolean isSupportConnected() {
304 return true;
305 }
306 }
307