    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    }