WeibullDistribution.java

  1. /*
  2.  * Licensed to the Apache Software Foundation (ASF) under one or more
  3.  * contributor license agreements.  See the NOTICE file distributed with
  4.  * this work for additional information regarding copyright ownership.
  5.  * The ASF licenses this file to You under the Apache License, Version 2.0
  6.  * (the "License"); you may not use this file except in compliance with
  7.  * the License.  You may obtain a copy of the License at
  8.  *
  9.  *      http://www.apache.org/licenses/LICENSE-2.0
  10.  *
  11.  * Unless required by applicable law or agreed to in writing, software
  12.  * distributed under the License is distributed on an "AS IS" BASIS,
  13.  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14.  * See the License for the specific language governing permissions and
  15.  * limitations under the License.
  16.  */

  18. package org.apache.commons.statistics.distribution;

  20. import org.apache.commons.numbers.gamma.LogGamma;
  21. import org.apache.commons.rng.UniformRandomProvider;
  22. import org.apache.commons.rng.sampling.distribution.ZigguratSampler;

  24. /**
  25.  * Implementation of the Weibull distribution.
  26.  *
  27.  * <p>The probability density function of \( X \) is:
  28.  *
  29.  * <p>\[ f(x;k,\lambda) = \frac{k}{\lambda}\left(\frac{x}{\lambda}\right)^{k-1}e^{-(x/\lambda)^{k}} \]
  30.  *
  31.  * <p>for \( k &gt; 0 \) the shape,
  32.  * \( \lambda &gt; 0 \) the scale, and
  33.  * \( x \in (0, \infty) \).
  34.  *
  35.  * <p>Note the special cases:
  36.  * <ul>
  37.  * <li>\( k = 1 \) is the exponential distribution
  38.  * <li>\( k = 2 \) is the Rayleigh distribution with scale \( \sigma = \frac {\lambda}{\sqrt{2}} \)
  39.  * </ul>
  40.  *
  41.  * @see <a href="https://en.wikipedia.org/wiki/Weibull_distribution">Weibull distribution (Wikipedia)</a>
  42.  * @see <a href="https://mathworld.wolfram.com/WeibullDistribution.html">Weibull distribution (MathWorld)</a>
  43.  */
  44. public final class WeibullDistribution extends AbstractContinuousDistribution {
  45.     /** Support lower bound. */
  46.     private static final double SUPPORT_LO = 0;
  47.     /** Support upper bound. */
  48.     private static final double SUPPORT_HI = Double.POSITIVE_INFINITY;
  49.     /** The shape parameter. */
  50.     private final double shape;
  51.     /** The scale parameter. */
  52.     private final double scale;
  53.     /** shape / scale. */
  54.     private final double shapeOverScale;
  55.     /** log(shape / scale). */
  56.     private final double logShapeOverScale;

  58.     /**
  59.      * @param shape Shape parameter.
  60.      * @param scale Scale parameter.
  61.      */
  62.     private WeibullDistribution(double shape,
  63.                                 double scale) {
  64.         this.scale = scale;
  65.         this.shape = shape;
  66.         shapeOverScale = shape / scale;
  67.         logShapeOverScale = Math.log(shapeOverScale);
  68.     }

  70.     /**
  71.      * Creates a Weibull distribution.
  72.      *
  73.      * @param shape Shape parameter.
  74.      * @param scale Scale parameter.
  75.      * @return the distribution
  76.      * @throws IllegalArgumentException if {@code shape <= 0} or {@code scale <= 0}.
  77.      */
  78.     public static WeibullDistribution of(double shape,
  79.                                          double scale) {
  80.         if (shape <= 0) {
  81.             throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
  82.                                             shape);
  83.         }
  84.         if (scale <= 0) {
  85.             throw new DistributionException(DistributionException.NOT_STRICTLY_POSITIVE,
  86.                                             scale);
  87.         }
  88.         return new WeibullDistribution(shape, scale);
  89.     }

  91.     /**
  92.      * Gets the shape parameter of this distribution.
  93.      *
  94.      * @return the shape parameter.
  95.      */
  96.     public double getShape() {
  97.         return shape;
  98.     }

  100.     /**
  101.      * Gets the scale parameter of this distribution.
  102.      *
  103.      * @return the scale parameter.
  104.      */
  105.     public double getScale() {
  106.         return scale;
  107.     }

  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.         }

  131.         final double xscale = x / scale;
  132.         final double xscalepow = Math.pow(xscale, shape - 1);

  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;

  141.         return shapeOverScale * xscalepow * Math.exp(-xscalepowshape);
  142.     }

  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.         }

  166.         final double xscale = x / scale;
  167.         final double logxscalepow = Math.log(xscale) * (shape - 1);

  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;

  177.         return logShapeOverScale + logxscalepow - xscalepowshape;
  178.     }

  180.     /** {@inheritDoc} */
  181.     @Override
  182.     public double cumulativeProbability(double x) {
  183.         if (x <= SUPPORT_LO) {
  184.             return 0;
  185.         }

  187.         return -Math.expm1(-Math.pow(x / scale, shape));
  188.     }

  190.     /** {@inheritDoc} */
  191.     @Override
  192.     public double survivalProbability(double x) {
  193.         if (x <= SUPPORT_LO) {
  194.             return 1;
  195.         }

  197.         return Math.exp(-Math.pow(x / scale, shape));
  198.     }

  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.     }

  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.     }

  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();

  248.         // Special case of exponential when shape is 1
  249.         return sh == 1 ? sc : sc * Math.exp(LogGamma.value(1 + (1 / sh)));
  250.     }

  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();

  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.     }

  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.     }

  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.     }

  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. }
Created with JaCoCo 0.8.12.202403310830Code Coverage Report for Apache Commons Statistics 1.1