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.MathRuntimeException;
023    import org.apache.commons.math.exception.util.LocalizedFormats;
024    import org.apache.commons.math.special.Gamma;
025    import org.apache.commons.math.util.FastMath;
026    
027    /**
028     * Default implementation of
029     * {@link org.apache.commons.math.distribution.WeibullDistribution}.
030     *
031     * @since 1.1
032     * @version $Revision: 1054524 $ $Date: 2011-01-03 05:59:18 +0100 (lun. 03 janv. 2011) $
033     */
034    public class WeibullDistributionImpl extends AbstractContinuousDistribution
035            implements WeibullDistribution, 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        /** The shape parameter. */
047        private double shape;
048    
049        /** The scale parameter. */
050        private double scale;
051    
052        /** Inverse cumulative probability accuracy */
053        private final double solverAbsoluteAccuracy;
054    
055        /** Cached numerical mean */
056        private double numericalMean = Double.NaN;
057    
058        /** Whether or not the numerical mean has been calculated */
059        private boolean numericalMeanIsCalculated = false;
060    
061        /** Cached numerical variance */
062        private double numericalVariance = Double.NaN;
063    
064        /** Whether or not the numerical variance has been calculated */
065        private boolean numericalVarianceIsCalculated = false;
066    
067        /**
068         * Creates weibull distribution with the given shape and scale and a
069         * location equal to zero.
070         * @param alpha the shape parameter.
071         * @param beta the scale parameter.
072         */
073        public WeibullDistributionImpl(double alpha, double beta){
074            this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
075        }
076    
077        /**
078         * Creates weibull distribution with the given shape, scale and inverse
079         * cumulative probability accuracy and a location equal to zero.
080         * @param alpha the shape parameter.
081         * @param beta the scale parameter.
082         * @param inverseCumAccuracy the maximum absolute error in inverse cumulative probability estimates
083         * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY})
084         * @since 2.1
085         */
086        public WeibullDistributionImpl(double alpha, double beta, double inverseCumAccuracy){
087            super();
088            setShapeInternal(alpha);
089            setScaleInternal(beta);
090            solverAbsoluteAccuracy = inverseCumAccuracy;
091        }
092    
093        /**
094         * For this distribution, X, this method returns P(X &lt; <code>x</code>).
095         * @param x the value at which the CDF is evaluated.
096         * @return CDF evaluated at <code>x</code>.
097         */
098        public double cumulativeProbability(double x) {
099            double ret;
100            if (x <= 0.0) {
101                ret = 0.0;
102            } else {
103                ret = 1.0 - FastMath.exp(-FastMath.pow(x / scale, shape));
104            }
105            return ret;
106        }
107    
108        /**
109         * Access the shape parameter.
110         * @return the shape parameter.
111         */
112        public double getShape() {
113            return shape;
114        }
115    
116        /**
117         * Access the scale parameter.
118         * @return the scale parameter.
119         */
120        public double getScale() {
121            return scale;
122        }
123    
124        /**
125         * Returns the probability density for a particular point.
126         *
127         * @param x The point at which the density should be computed.
128         * @return The pdf at point x.
129         * @since 2.1
130         */
131        @Override
132        public double density(double x) {
133            if (x < 0) {
134                return 0;
135            }
136    
137            final double xscale = x / scale;
138            final double xscalepow = FastMath.pow(xscale, shape - 1);
139    
140            /*
141             * FastMath.pow(x / scale, shape) =
142             * FastMath.pow(xscale, shape) =
143             * FastMath.pow(xscale, shape - 1) * xscale
144             */
145            final double xscalepowshape = xscalepow * xscale;
146    
147            return (shape / scale) * xscalepow * FastMath.exp(-xscalepowshape);
148        }
149    
150        /**
151         * For this distribution, X, this method returns the critical point x, such
152         * that P(X &lt; x) = <code>p</code>.
153         * <p>
154         * Returns <code>Double.NEGATIVE_INFINITY</code> for p=0 and
155         * <code>Double.POSITIVE_INFINITY</code> for p=1.</p>
156         *
157         * @param p the desired probability
158         * @return x, such that P(X &lt; x) = <code>p</code>
159         * @throws IllegalArgumentException if <code>p</code> is not a valid
160         *         probability.
161         */
162        @Override
163        public double inverseCumulativeProbability(double p) {
164            double ret;
165            if (p < 0.0 || p > 1.0) {
166                throw MathRuntimeException.createIllegalArgumentException(
167                      LocalizedFormats.OUT_OF_RANGE_SIMPLE, p, 0.0, 1.0);
168            } else if (p == 0) {
169                ret = 0.0;
170            } else  if (p == 1) {
171                ret = Double.POSITIVE_INFINITY;
172            } else {
173                ret = scale * FastMath.pow(-FastMath.log(1.0 - p), 1.0 / shape);
174            }
175            return ret;
176        }
177    
178        /**
179         * Modify the shape parameter.
180         * @param alpha the new shape parameter value.
181         * @deprecated as of 2.1 (class will become immutable in 3.0)
182         */
183        @Deprecated
184        public void setShape(double alpha) {
185            setShapeInternal(alpha);
186            invalidateParameterDependentMoments();
187        }
188        /**
189         * Modify the shape parameter.
190         * @param alpha the new shape parameter value.
191         */
192        private void setShapeInternal(double alpha) {
193            if (alpha <= 0.0) {
194                throw MathRuntimeException.createIllegalArgumentException(
195                      LocalizedFormats.NOT_POSITIVE_SHAPE,
196                      alpha);
197            }
198            this.shape = alpha;
199        }
200    
201        /**
202         * Modify the scale parameter.
203         * @param beta the new scale parameter value.
204         * @deprecated as of 2.1 (class will become immutable in 3.0)
205         */
206        @Deprecated
207        public void setScale(double beta) {
208            setScaleInternal(beta);
209            invalidateParameterDependentMoments();
210        }
211        /**
212         * Modify the scale parameter.
213         * @param beta the new scale parameter value.
214         */
215        private void setScaleInternal(double beta) {
216            if (beta <= 0.0) {
217                throw MathRuntimeException.createIllegalArgumentException(
218                      LocalizedFormats.NOT_POSITIVE_SCALE,
219                      beta);
220            }
221            this.scale = beta;
222        }
223    
224        /**
225         * Access the domain value lower bound, based on <code>p</code>, used to
226         * bracket a CDF root.  This method is used by
227         * {@link #inverseCumulativeProbability(double)} to find critical values.
228         *
229         * @param p the desired probability for the critical value
230         * @return domain value lower bound, i.e.
231         *         P(X &lt; <i>lower bound</i>) &lt; <code>p</code>
232         */
233        @Override
234        protected double getDomainLowerBound(double p) {
235            return 0.0;
236        }
237    
238        /**
239         * Access the domain value upper bound, based on <code>p</code>, used to
240         * bracket a CDF root.  This method is used by
241         * {@link #inverseCumulativeProbability(double)} to find critical values.
242         *
243         * @param p the desired probability for the critical value
244         * @return domain value upper bound, i.e.
245         *         P(X &lt; <i>upper bound</i>) &gt; <code>p</code>
246         */
247        @Override
248        protected double getDomainUpperBound(double p) {
249            return Double.MAX_VALUE;
250        }
251    
252        /**
253         * Access the initial domain value, based on <code>p</code>, used to
254         * bracket a CDF root.  This method is used by
255         * {@link #inverseCumulativeProbability(double)} to find critical values.
256         *
257         * @param p the desired probability for the critical value
258         * @return initial domain value
259         */
260        @Override
261        protected double getInitialDomain(double p) {
262            // use median
263            return FastMath.pow(scale * FastMath.log(2.0), 1.0 / shape);
264        }
265    
266        /**
267         * Return the absolute accuracy setting of the solver used to estimate
268         * inverse cumulative probabilities.
269         *
270         * @return the solver absolute accuracy
271         * @since 2.1
272         */
273        @Override
274        protected double getSolverAbsoluteAccuracy() {
275            return solverAbsoluteAccuracy;
276        }
277    
278        /**
279         * Returns the lower bound of the support for the distribution.
280         *
281         * The lower bound of the support is always 0 no matter the parameters.
282         *
283         * @return lower bound of the support (always 0)
284         * @since 2.2
285         */
286        public double getSupportLowerBound() {
287            return 0;
288        }
289    
290        /**
291         * Returns the upper bound of the support for the distribution.
292         *
293         * The upper bound of the support is always positive infinity
294         * no matter the parameters.
295         *
296         * @return upper bound of the support (always Double.POSITIVE_INFINITY)
297         * @since 2.2
298         */
299        public double getSupportUpperBound() {
300            return Double.POSITIVE_INFINITY;
301        }
302    
303        /**
304         * Calculates the mean.
305         *
306         * The mean is <code>scale * Gamma(1 + (1 / shape))</code>
307         * where <code>Gamma(...)</code> is the Gamma-function
308         *
309         * @return the mean
310         * @since 2.2
311         */
312        protected double calculateNumericalMean() {
313            final double sh = getShape();
314            final double sc = getScale();
315    
316            return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
317        }
318    
319        /**
320         * Calculates the variance.
321         *
322         * The variance is
323         * <code>scale^2 * Gamma(1 + (2 / shape)) - mean^2</code>
324         * where <code>Gamma(...)</code> is the Gamma-function
325         *
326         * @return the variance
327         * @since 2.2
328         */
329        private double calculateNumericalVariance() {
330            final double sh = getShape();
331            final double sc = getScale();
332            final double mn = getNumericalMean();
333    
334            return (sc * sc) *
335                FastMath.exp(Gamma.logGamma(1 + (2 / sh))) -
336                (mn * mn);
337        }
338    
339        /**
340         * Returns the mean of the distribution.
341         *
342         * @return the mean or Double.NaN if it's not defined
343         * @since 2.2
344         */
345        public double getNumericalMean() {
346            if (!numericalMeanIsCalculated) {
347                numericalMean = calculateNumericalMean();
348                numericalMeanIsCalculated = true;
349            }
350    
351            return numericalMean;
352        }
353    
354        /**
355         * Returns the variance of the distribution.
356         *
357         * @return the variance (possibly Double.POSITIVE_INFINITY as
358         * for certain cases in {@link TDistributionImpl}) or
359         * Double.NaN if it's not defined
360         * @since 2.2
361         */
362        public double getNumericalVariance() {
363            if (!numericalVarianceIsCalculated) {
364                numericalVariance = calculateNumericalVariance();
365                numericalVarianceIsCalculated = true;
366            }
367    
368            return numericalVariance;
369        }
370    
371        /**
372         * Invalidates the cached mean and variance.
373         */
374        private void invalidateParameterDependentMoments() {
375            numericalMeanIsCalculated = false;
376            numericalVarianceIsCalculated = false;
377        }
378    }