Source code

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;
021import org.apache.commons.rng.UniformRandomProvider;
022import org.apache.commons.rng.sampling.distribution.ZigguratSampler;
023
024/**
025 * Implementation of the Weibull distribution.
026 *
027 * <p>The probability density function of \( X \) is:
028 *
029 * <p>\[ f(x;k,\lambda) = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1}e^{-(x/\lambda)^{k}} \]
030 *
031 * <p>for \( k &gt; 0 \) the shape,
032 * \( \lambda &gt; 0 \) the scale, and
033 * \( x \in (0, \infty) \).
034 *
035 * <p>Note the special cases:
036 * <ul>
037 * <li>\( k = 1 \) is the exponential distribution
038 * <li>\( k = 2 \) is the Rayleigh distribution with scale \( \sigma = \frac {\lambda}{\sqrt{2}} \)
039 * </ul>
040 *
041 * @see <a href="https://en.wikipedia.org/wiki/Weibull_distribution">Weibull distribution (Wikipedia)</a>
042 * @see <a href="https://mathworld.wolfram.com/WeibullDistribution.html">Weibull distribution (MathWorld)</a>
043 */
044public final class WeibullDistribution extends AbstractContinuousDistribution {
045    /** Support lower bound. */
046    private static final double SUPPORT_LO = 0;
047    /** Support upper bound. */
048    private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
049    /** The shape parameter. */
050    private final double shape;
051    /** The scale parameter. */
052    private final double scale;
053    /** shape / scale. */
054    private final double shapeOverScale;
055    /** log(shape / scale). */
056    private final double logShapeOverScale;
057
058    /**
059     * @param shape Shape parameter.
060     * @param scale Scale parameter.
061     */
062    private WeibullDistribution(double shape,
063                                double scale) {
064        this.scale = scale;
065        this.shape = shape;
066        shapeOverScale = shape / scale;
067        logShapeOverScale = Math.log(shapeOverScale);
068    }
069
070    /**
071     * Creates a Weibull distribution.
072     *
073     * @param shape Shape parameter.
074     * @param scale Scale parameter.
075     * @return the distribution
076     * @throws IllegalArgumentException if {@code shape <= 0} or {@code scale <= 0}.
077     */
078    public static WeibullDistribution of(double shape,
079                                         double scale) {
080        if (shape <= 0) {
081            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
082                                            shape);
083        }
084        if (scale <= 0) {
085            throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
086                                            scale);
087        }
088        return new WeibullDistribution(shape, scale);
089    }
090
091    /**
092     * Gets the shape parameter of this distribution.
093     *
094     * @return the shape parameter.
095     */
096    public double getShape() {
097        return shape;
098    }
099
100    /**
101     * Gets the scale parameter of this distribution.
102     *
103     * @return the scale parameter.
104     */
105    public double getScale() {
106        return scale;
107    }
108
109    /** {@inheritDoc}
110     *
111     * <p>Returns the limit when {@code x = 0}:
112     * <ul>
113     * <li>{@code shape < 1}: Infinity
114     * <li>{@code shape == 1}: 1 / scale
115     * <li>{@code shape > 1}: 0
116     * </ul>
117     */
118    @Override
119    public double density(double x) {
120        if (x <= SUPPORT_LO || x >= SUPPORT_HI) {
121            // Special case x=0
122            if (x == SUPPORT_LO && shape <= 1) {
123                return shape == 1 ?
124                    // Exponential distribution
125                    shapeOverScale :
126                    Double.POSITIVE_INFINITY;
127            }
128            return 0;
129        }
130
131        final double xscale = x / scale;
132        final double xscalepow = Math.pow(xscale, shape - 1);
133
134        /*
135         * Math.pow(x / scale, shape) =
136         * Math.pow(xscale, shape) =
137         * Math.pow(xscale, shape - 1) * xscale
138         */
139        final double xscalepowshape = xscalepow * xscale;
140
141        return shapeOverScale * xscalepow * Math.exp(-xscalepowshape);
142    }
143
144    /** {@inheritDoc}
145     *
146     * <p>Returns the limit when {@code x = 0}:
147     * <ul>
148     * <li>{@code shape < 1}: Infinity
149     * <li>{@code shape == 1}: log(1 / scale)
150     * <li>{@code shape > 1}: -Infinity
151     * </ul>
152     */
153    @Override
154    public double logDensity(double x) {
155        if (x <= SUPPORT_LO || x >= SUPPORT_HI) {
156            // Special case x=0
157            if (x == SUPPORT_LO && shape <= 1) {
158                return shape == 1 ?
159                    // Exponential distribution
160                    logShapeOverScale :
161                    Double.POSITIVE_INFINITY;
162            }
163            return Double.NEGATIVE_INFINITY;
164        }
165
166        final double xscale = x / scale;
167        final double logxscalepow = Math.log(xscale) * (shape - 1);
168
169        /*
170         * Math.pow(x / scale, shape) =
171         * Math.pow(xscale, shape) =
172         * Math.pow(xscale, shape - 1) * xscale
173         * Math.exp(log(xscale) * (shape - 1)) * xscale
174         */
175        final double xscalepowshape = Math.exp(logxscalepow) * xscale;
176
177        return logShapeOverScale + logxscalepow - xscalepowshape;
178    }
179
180    /** {@inheritDoc} */
181    @Override
182    public double cumulativeProbability(double x) {
183        if (x <= SUPPORT_LO) {
184            return 0;
185        }
186
187        return -Math.expm1(-Math.pow(x / scale, shape));
188    }
189
190    /** {@inheritDoc} */
191    @Override
192    public double survivalProbability(double x) {
193        if (x <= SUPPORT_LO) {
194            return 1;
195        }
196
197        return Math.exp(-Math.pow(x / scale, shape));
198    }
199
200    /**
201     * {@inheritDoc}
202     *
203     * <p>Returns {@code 0} when {@code p == 0} and
204     * {@link Double#POSITIVE_INFINITY} when {@code p == 1}.
205     */
206    @Override
207    public double inverseCumulativeProbability(double p) {
208        ArgumentUtils.checkProbability(p);
209        if (p == 0) {
210            return 0.0;
211        } else  if (p == 1) {
212            return Double.POSITIVE_INFINITY;
213        }
214        return scale * Math.pow(-Math.log1p(-p), 1.0 / shape);
215    }
216
217    /**
218     * {@inheritDoc}
219     *
220     * <p>Returns {@code 0} when {@code p == 1} and
221     * {@link Double#POSITIVE_INFINITY} when {@code p == 0}.
222     */
223    @Override
224    public double inverseSurvivalProbability(double p) {
225        ArgumentUtils.checkProbability(p);
226        if (p == 1) {
227            return 0.0;
228        } else  if (p == 0) {
229            return Double.POSITIVE_INFINITY;
230        }
231        return scale * Math.pow(-Math.log(p), 1.0 / shape);
232    }
233
234    /**
235     * {@inheritDoc}
236     *
237     * <p>For shape parameter \( k \) and scale parameter \( \lambda \), the mean is:
238     *
239     * <p>\[ \lambda \, \Gamma(1+\frac{1}{k}) \]
240     *
241     * <p>where \( \Gamma \) is the Gamma-function.
242     */
243    @Override
244    public double getMean() {
245        final double sh = getShape();
246        final double sc = getScale();
247
248        // Special case of exponential when shape is 1
249        return sh == 1 ? sc : sc * Math.exp(LogGamma.value(1 + (1 / sh)));
250    }
251
252    /**
253     * {@inheritDoc}
254     *
255     * <p>For shape parameter \( k \) and scale parameter \( \lambda \), the variance is:
256     *
257     * <p>\[ \lambda^2 \left[ \Gamma\left(1+\frac{2}{k}\right) -
258     *                        \left(\Gamma\left(1+\frac{1}{k}\right)\right)^2 \right] \]
259     *
260     * <p>where \( \Gamma \) is the Gamma-function.
261     */
262    @Override
263    public double getVariance() {
264        final double sh = getShape();
265        final double sc = getScale();
266        final double mn = getMean();
267
268        // Special case of exponential when shape is 1
269        return sh == 1 ?
270               sc * sc :
271               (sc * sc) * Math.exp(LogGamma.value(1 + (2 / sh))) -
272               (mn * mn);
273    }
274
275    /**
276     * {@inheritDoc}
277     *
278     * <p>The lower bound of the support is always 0.
279     *
280     * @return 0.
281     */
282    @Override
283    public double getSupportLowerBound() {
284        return SUPPORT_LO;
285    }
286
287    /**
288     * {@inheritDoc}
289     *
290     * <p>The upper bound of the support is always positive infinity.
291     *
292     * @return {@linkplain Double#POSITIVE_INFINITY positive infinity}.
293     */
294    @Override
295    public double getSupportUpperBound() {
296        return SUPPORT_HI;
297    }
298
299    /** {@inheritDoc} */
300    @Override
301    public ContinuousDistribution.Sampler createSampler(final UniformRandomProvider rng) {
302        // Special case: shape=1 is the exponential distribution
303        if (shape == 1) {
304            // Exponential distribution sampler.
305            return ZigguratSampler.Exponential.of(rng, scale)::sample;
306        }
307        return super.createSampler(rng);
308    }
309}