    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