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
018package org.apache.commons.statistics.distribution;
019
020import org.apache.commons.numbers.gamma.LogGamma;
021
022/**
023 * Implementation of the Weibull distribution. This implementation uses the
024 * two parameter form of the distribution defined by
025 * <a href="http://mathworld.wolfram.com/WeibullDistribution.html">
026 * Weibull Distribution</a>, equations (1) and (2).
027 *
028 * @see <a href="http://en.wikipedia.org/wiki/Weibull_distribution">Weibull distribution (Wikipedia)</a>
029 * @see <a href="http://mathworld.wolfram.com/WeibullDistribution.html">Weibull distribution (MathWorld)</a>
030 */
031public class WeibullDistribution extends AbstractContinuousDistribution {
032    /** Support lower bound. */
033    private static final double SUPPORT_LO = 0;
034    /** Support upper bound. */
035    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
036    /** The shape parameter. */
037    private final double shape;
038    /** The scale parameter. */
039    private final double scale;
040
041    /**
042     * Creates a distribution.
043     *
044     * @param alpha Shape parameter.
045     * @param beta Scale parameter.
046     * @throws IllegalArgumentException if {@code alpha <= 0} or {@code beta <= 0}.
047     */
048    public WeibullDistribution(double alpha,
049                               double beta) {
050        if (alpha <= 0) {
051            throw new DistributionException(DistributionException.NEGATIVE,
052                                            alpha);
053        }
054        if (beta <= 0) {
055            throw new DistributionException(DistributionException.NEGATIVE,
056                                            beta);
057        }
058        scale = beta;
059        shape = alpha;
060    }
061
062    /**
063     * Access the shape parameter, {@code alpha}.
064     *
065     * @return the shape parameter, {@code alpha}.
066     */
067    public double getShape() {
068        return shape;
069    }
070
071    /**
072     * Access the scale parameter, {@code beta}.
073     *
074     * @return the scale parameter, {@code beta}.
075     */
076    public double getScale() {
077        return scale;
078    }
079
080    /** {@inheritDoc} */
081    @Override
082    public double density(double x) {
083        if (x <= SUPPORT_LO ||
084            x >= SUPPORT_HI) {
085            return 0;
086        }
087
088        final double xscale = x / scale;
089        final double xscalepow = Math.pow(xscale, shape - 1);
090
091        /*
092         * Math.pow(x / scale, shape) =
093         * Math.pow(xscale, shape) =
094         * Math.pow(xscale, shape - 1) * xscale
095         */
096        final double xscalepowshape = xscalepow * xscale;
097
098        return (shape / scale) * xscalepow * Math.exp(-xscalepowshape);
099    }
100
101    /** {@inheritDoc} */
102    @Override
103    public double logDensity(double x) {
104        if (x <= SUPPORT_LO ||
105            x >= SUPPORT_HI) {
106            return Double.NEGATIVE_INFINITY;
107        }
108
109        final double xscale = x / scale;
110        final double logxscalepow = Math.log(xscale) * (shape - 1);
111
112        /*
113         * Math.pow(x / scale, shape) =
114         * Math.pow(xscale, shape) =
115         * Math.pow(xscale, shape - 1) * xscale
116         */
117        final double xscalepowshape = Math.exp(logxscalepow) * xscale;
118
119        return Math.log(shape / scale) + logxscalepow - xscalepowshape;
120    }
121
122    /** {@inheritDoc} */
123    @Override
124    public double cumulativeProbability(double x) {
125        if (x <= SUPPORT_LO) {
126            return 0;
127        }
128
129        return 1 - Math.exp(-Math.pow(x / scale, shape));
130    }
131
132    /**
133     * {@inheritDoc}
134     *
135     * Returns {@code 0} when {@code p == 0} and
136     * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
137     */
138    @Override
139    public double inverseCumulativeProbability(double p) {
140        double ret;
141        if (p < 0 ||
142            p > 1) {
143            throw new DistributionException(DistributionException.INVALID_PROBABILITY, p);
144        } else if (p == 0) {
145            ret = 0.0;
146        } else  if (p == 1) {
147            ret = Double.POSITIVE_INFINITY;
148        } else {
149            ret = scale * Math.pow(-Math.log1p(-p), 1.0 / shape);
150        }
151        return ret;
152    }
153
154    /**
155     * {@inheritDoc}
156     *
157     * The mean is {@code scale * Gamma(1 + (1 / shape))}, where {@code Gamma()}
158     * is the Gamma-function.
159     */
160    @Override
161    public double getMean() {
162        final double sh = getShape();
163        final double sc = getScale();
164
165        return sc * Math.exp(LogGamma.value(1 + (1 / sh)));
166    }
167
168    /**
169     * {@inheritDoc}
170     *
171     * The variance is {@code scale^2 * Gamma(1 + (2 / shape)) - mean^2}
172     * where {@code Gamma()} is the Gamma-function.
173     */
174    @Override
175    public double getVariance() {
176        final double sh = getShape();
177        final double sc = getScale();
178        final double mn = getMean();
179
180        return (sc * sc) * Math.exp(LogGamma.value(1 + (2 / sh))) -
181               (mn * mn);
182    }
183
184    /**
185     * {@inheritDoc}
186     *
187     * The lower bound of the support is always 0 no matter the parameters.
188     *
189     * @return lower bound of the support (always 0)
190     */
191    @Override
192    public double getSupportLowerBound() {
193        return SUPPORT_LO;
194    }
195
196    /**
197     * {@inheritDoc}
198     *
199     * The upper bound of the support is always positive infinity
200     * no matter the parameters.
201     *
202     * @return upper bound of the support (always
203     * {@code Double.POSITIVE_INFINITY})
204     */
205    @Override
206    public double getSupportUpperBound() {
207        return SUPPORT_HI;
208    }
209
210    /**
211     * {@inheritDoc}
212     *
213     * The support of this distribution is connected.
214     *
215     * @return {@code true}
216     */
217    @Override
218    public boolean isSupportConnected() {
219        return true;
220    }
221}
222