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.exception.OutOfRangeException;
023 import org.apache.commons.math.exception.NotStrictlyPositiveException;
024 import org.apache.commons.math.exception.util.LocalizedFormats;
025 import org.apache.commons.math.special.Gamma;
026 import org.apache.commons.math.util.FastMath;
027
028 /**
029 * Default implementation of
030 * {@link org.apache.commons.math.distribution.WeibullDistribution}.
031 *
032 * @since 1.1
033 * @version $Id: WeibullDistributionImpl.java 1131229 2011-06-03 20:49:25Z luc $
034 */
035 public class WeibullDistributionImpl extends AbstractContinuousDistribution
036 implements WeibullDistribution, Serializable {
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 /** The shape parameter. */
045 private final double shape;
046 /** The scale parameter. */
047 private final double scale;
048 /** Inverse cumulative probability accuracy. */
049 private final double solverAbsoluteAccuracy;
050
051 /**
052 * Create a Weibull distribution with the given shape and scale and a
053 * location equal to zero.
054 *
055 * @param alpha Shape parameter.
056 * @param beta Scale parameter.
057 */
058 public WeibullDistributionImpl(double alpha, double beta) {
059 this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
060 }
061
062 /**
063 * Create a Weibull distribution with the given shape, scale and inverse
064 * cumulative probability accuracy and a location equal to zero.
065 *
066 * @param alpha Shape parameter.
067 * @param beta Scale parameter.
068 * @param inverseCumAccuracy Maximum absolute error in inverse
069 * cumulative probability estimates
070 * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
071 * @throws NotStrictlyPositiveException if {@code alpha <= 0} or
072 * {@code beta <= 0}.
073 * @since 2.1
074 */
075 public WeibullDistributionImpl(double alpha, double beta,
076 double inverseCumAccuracy) {
077 if (alpha <= 0) {
078 throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE,
079 alpha);
080 }
081 if (beta <= 0) {
082 throw new NotStrictlyPositiveException(LocalizedFormats.SCALE,
083 beta);
084 }
085 scale = beta;
086 shape = alpha;
087 solverAbsoluteAccuracy = inverseCumAccuracy;
088 }
089
090 /**
091 * For this distribution, {@code X}, this method returns {@code P(X < x)}.
092 *
093 * @param x Value at which the CDF is evaluated.
094 * @return the CDF evaluated at {@code x}.
095 */
096 public double cumulativeProbability(double x) {
097 double ret;
098 if (x <= 0.0) {
099 ret = 0.0;
100 } else {
101 ret = 1.0 - FastMath.exp(-FastMath.pow(x / scale, shape));
102 }
103 return ret;
104 }
105
106 /**
107 * {@inheritDoc}
108 */
109 public double getShape() {
110 return shape;
111 }
112
113 /**
114 * {@inheritDoc}
115 */
116 public double getScale() {
117 return scale;
118 }
119
120 /**
121 * {@inheritDoc}
122 */
123 @Override
124 public double density(double x) {
125 if (x < 0) {
126 return 0;
127 }
128
129 final double xscale = x / scale;
130 final double xscalepow = FastMath.pow(xscale, shape - 1);
131
132 /*
133 * FastMath.pow(x / scale, shape) =
134 * FastMath.pow(xscale, shape) =
135 * FastMath.pow(xscale, shape - 1) * xscale
136 */
137 final double xscalepowshape = xscalepow * xscale;
138
139 return (shape / scale) * xscalepow * FastMath.exp(-xscalepowshape);
140 }
141
142 /**
143 * For this distribution, {@code X}, this method returns the critical
144 * point {@code x}, such that {@code P(X < x) = p}.
145 * It will return {@code Double.NEGATIVE_INFINITY} when p = 0 and
146 * {@code Double.POSITIVE_INFINITY} when p = 1.
147 *
148 * @param p Desired probability.
149 * @return {@code x}, such that {@code P(X < x) = p}.
150 * @throws OutOfRangeException if {@code p} is not a valid probability.
151 */
152 @Override
153 public double inverseCumulativeProbability(double p) {
154 double ret;
155 if (p < 0.0 || p > 1.0) {
156 throw new OutOfRangeException(p, 0.0, 1.0);
157 } else if (p == 0) {
158 ret = 0.0;
159 } else if (p == 1) {
160 ret = Double.POSITIVE_INFINITY;
161 } else {
162 ret = scale * FastMath.pow(-FastMath.log(1.0 - p), 1.0 / shape);
163 }
164 return ret;
165 }
166
167
168 /**
169 * Access the domain value lower bound, based on {@code p}, used to
170 * bracket a CDF root. This method is used by
171 * {@link #inverseCumulativeProbability(double)} to find critical values.
172 *
173 * @param p Desired probability for the critical value.
174 * @return the domain value lower bound, i.e. {@code P(X < 'lower bound') < p}.
175 */
176 @Override
177 protected double getDomainLowerBound(double p) {
178 return 0;
179 }
180
181 /**
182 * Access the domain value upper bound, based on {@code p}, used to
183 * bracket a CDF root. This method is used by
184 * {@link #inverseCumulativeProbability(double)} to find critical values.
185 *
186 * @param p Desired probability for the critical value.
187 * @return the domain value upper bound, i.e. {@code P(X < 'upper bound') > p}.
188 */
189 @Override
190 protected double getDomainUpperBound(double p) {
191 return Double.MAX_VALUE;
192 }
193
194 /**
195 * Access the initial domain value, based on {@code p}, used to
196 * bracket a CDF root. This method is used by
197 * {@link #inverseCumulativeProbability(double)} to find critical values.
198 *
199 * @param p Desired probability for the critical value.
200 * @return the initial domain value.
201 */
202 @Override
203 protected double getInitialDomain(double p) {
204 // use median
205 return FastMath.pow(scale * FastMath.log(2.0), 1.0 / shape);
206 }
207
208 /**
209 * Return the absolute accuracy setting of the solver used to estimate
210 * inverse cumulative probabilities.
211 *
212 * @return the solver absolute accuracy.
213 * @since 2.1
214 */
215 @Override
216 protected double getSolverAbsoluteAccuracy() {
217 return solverAbsoluteAccuracy;
218 }
219
220 /**
221 * {@inheritDoc}
222 *
223 * The lower bound of the support is always 0 no matter the parameters.
224 *
225 * @return lower bound of the support (always 0)
226 */
227 @Override
228 public double getSupportLowerBound() {
229 return 0;
230 }
231
232 /**
233 * {@inheritDoc}
234 *
235 * The upper bound of the support is always positive infinity
236 * no matter the parameters.
237 *
238 * @return upper bound of the support (always Double.POSITIVE_INFINITY)
239 */
240 @Override
241 public double getSupportUpperBound() {
242 return Double.POSITIVE_INFINITY;
243 }
244
245 /**
246 * {@inheritDoc}
247 *
248 * The mean is <code>scale * Gamma(1 + (1 / shape))</code>
249 * where <code>Gamma(...)</code> is the Gamma-function
250 *
251 * @return {@inheritDoc}
252 */
253 @Override
254 protected double calculateNumericalMean() {
255 final double sh = getShape();
256 final double sc = getScale();
257
258 return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
259 }
260
261 /**
262 * {@inheritDoc}
263 *
264 * The variance is
265 * <code>scale^2 * Gamma(1 + (2 / shape)) - mean^2</code>
266 * where <code>Gamma(...)</code> is the Gamma-function
267 *
268 * @return {@inheritDoc}
269 */
270 @Override
271 protected double calculateNumericalVariance() {
272 final double sh = getShape();
273 final double sc = getScale();
274 final double mn = getNumericalMean();
275
276 return (sc * sc) *
277 FastMath.exp(Gamma.logGamma(1 + (2 / sh))) -
278 (mn * mn);
279 }
280
281 /**
282 * {@inheritDoc}
283 */
284 @Override
285 public boolean isSupportLowerBoundInclusive() {
286 return true;
287 }
288
289 /**
290 * {@inheritDoc}
291 */
292 @Override
293 public boolean isSupportUpperBoundInclusive() {
294 return false;
295 }
296 }