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.math3.distribution;

  20. import org.apache.commons.math3.exception.NotStrictlyPositiveException;
  21. import org.apache.commons.math3.exception.OutOfRangeException;
  22. import org.apache.commons.math3.exception.util.LocalizedFormats;
  23. import org.apache.commons.math3.random.RandomGenerator;
  24. import org.apache.commons.math3.random.Well19937c;
  25. import org.apache.commons.math3.special.Gamma;
  26. import org.apache.commons.math3.util.FastMath;

  28. /**
  29.  * Implementation of the Weibull distribution. This implementation uses the
  30.  * two parameter form of the distribution defined by
  31.  * <a href="http://mathworld.wolfram.com/WeibullDistribution.html">
  32.  * Weibull Distribution</a>, equations (1) and (2).
  33.  *
  34.  * @see <a href="http://en.wikipedia.org/wiki/Weibull_distribution">Weibull distribution (Wikipedia)</a>
  35.  * @see <a href="http://mathworld.wolfram.com/WeibullDistribution.html">Weibull distribution (MathWorld)</a>
  36.  * @since 1.1 (changed to concrete class in 3.0)
  37.  */
  38. public class WeibullDistribution extends AbstractRealDistribution {
  39.     /**
  40.      * Default inverse cumulative probability accuracy.
  41.      * @since 2.1
  42.      */
  43.     public static final double DEFAULT_INVERSE_ABSOLUTE_ACCURACY = 1e-9;
  44.     /** Serializable version identifier. */
  45.     private static final long serialVersionUID = 8589540077390120676L;
  46.     /** The shape parameter. */
  47.     private final double shape;
  48.     /** The scale parameter. */
  49.     private final double scale;
  50.     /** Inverse cumulative probability accuracy. */
  51.     private final double solverAbsoluteAccuracy;
  52.     /** Cached numerical mean */
  53.     private double numericalMean = Double.NaN;
  54.     /** Whether or not the numerical mean has been calculated */
  55.     private boolean numericalMeanIsCalculated = false;
  56.     /** Cached numerical variance */
  57.     private double numericalVariance = Double.NaN;
  58.     /** Whether or not the numerical variance has been calculated */
  59.     private boolean numericalVarianceIsCalculated = false;

  61.     /**
  62.      * Create a Weibull distribution with the given shape and scale and a
  63.      * location equal to zero.
  64.      * <p>
  65.      * <b>Note:</b> this constructor will implicitly create an instance of
  66.      * {@link Well19937c} as random generator to be used for sampling only (see
  67.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  68.      * needed for the created distribution, it is advised to pass {@code null}
  69.      * as random generator via the appropriate constructors to avoid the
  70.      * additional initialisation overhead.
  71.      *
  72.      * @param alpha Shape parameter.
  73.      * @param beta Scale parameter.
  74.      * @throws NotStrictlyPositiveException if {@code alpha <= 0} or
  75.      * {@code beta <= 0}.
  76.      */
  77.     public WeibullDistribution(double alpha, double beta)
  78.         throws NotStrictlyPositiveException {
  79.         this(alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  80.     }

  82.     /**
  83.      * Create a Weibull distribution with the given shape, scale and inverse
  84.      * cumulative probability accuracy and a location equal to zero.
  85.      * <p>
  86.      * <b>Note:</b> this constructor will implicitly create an instance of
  87.      * {@link Well19937c} as random generator to be used for sampling only (see
  88.      * {@link #sample()} and {@link #sample(int)}). In case no sampling is
  89.      * needed for the created distribution, it is advised to pass {@code null}
  90.      * as random generator via the appropriate constructors to avoid the
  91.      * additional initialisation overhead.
  92.      *
  93.      * @param alpha Shape parameter.
  94.      * @param beta Scale parameter.
  95.      * @param inverseCumAccuracy Maximum absolute error in inverse
  96.      * cumulative probability estimates
  97.      * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
  98.      * @throws NotStrictlyPositiveException if {@code alpha <= 0} or
  99.      * {@code beta <= 0}.
  100.      * @since 2.1
  101.      */
  102.     public WeibullDistribution(double alpha, double beta,
  103.                                double inverseCumAccuracy) {
  104.         this(new Well19937c(), alpha, beta, inverseCumAccuracy);
  105.     }

  107.     /**
  108.      * Creates a Weibull distribution.
  109.      *
  110.      * @param rng Random number generator.
  111.      * @param alpha Shape parameter.
  112.      * @param beta Scale parameter.
  113.      * @throws NotStrictlyPositiveException if {@code alpha <= 0} or {@code beta <= 0}.
  114.      * @since 3.3
  115.      */
  116.     public WeibullDistribution(RandomGenerator rng, double alpha, double beta)
  117.         throws NotStrictlyPositiveException {
  118.         this(rng, alpha, beta, DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
  119.     }

  121.     /**
  122.      * Creates a Weibull distribution.
  123.      *
  124.      * @param rng Random number generator.
  125.      * @param alpha Shape parameter.
  126.      * @param beta Scale parameter.
  127.      * @param inverseCumAccuracy Maximum absolute error in inverse
  128.      * cumulative probability estimates
  129.      * (defaults to {@link #DEFAULT_INVERSE_ABSOLUTE_ACCURACY}).
  130.      * @throws NotStrictlyPositiveException if {@code alpha <= 0} or {@code beta <= 0}.
  131.      * @since 3.1
  132.      */
  133.     public WeibullDistribution(RandomGenerator rng,
  134.                                double alpha,
  135.                                double beta,
  136.                                double inverseCumAccuracy)
  137.         throws NotStrictlyPositiveException {
  138.         super(rng);

  140.         if (alpha <= 0) {
  141.             throw new NotStrictlyPositiveException(LocalizedFormats.SHAPE,
  142.                                                    alpha);
  143.         }
  144.         if (beta <= 0) {
  145.             throw new NotStrictlyPositiveException(LocalizedFormats.SCALE,
  146.                                                    beta);
  147.         }
  148.         scale = beta;
  149.         shape = alpha;
  150.         solverAbsoluteAccuracy = inverseCumAccuracy;
  151.     }

  153.     /**
  154.      * Access the shape parameter, {@code alpha}.
  155.      *
  156.      * @return the shape parameter, {@code alpha}.
  157.      */
  158.     public double getShape() {
  159.         return shape;
  160.     }

  162.     /**
  163.      * Access the scale parameter, {@code beta}.
  164.      *
  165.      * @return the scale parameter, {@code beta}.
  166.      */
  167.     public double getScale() {
  168.         return scale;
  169.     }

  171.     /** {@inheritDoc} */
  172.     public double density(double x) {
  173.         if (x < 0) {
  174.             return 0;
  175.         }

  177.         final double xscale = x / scale;
  178.         final double xscalepow = FastMath.pow(xscale, shape - 1);

  180.         /*
  181.          * FastMath.pow(x / scale, shape) =
  182.          * FastMath.pow(xscale, shape) =
  183.          * FastMath.pow(xscale, shape - 1) * xscale
  184.          */
  185.         final double xscalepowshape = xscalepow * xscale;

  187.         return (shape / scale) * xscalepow * FastMath.exp(-xscalepowshape);
  188.     }

  190.     /** {@inheritDoc} */
  191.     @Override
  192.     public double logDensity(double x) {
  193.         if (x < 0) {
  194.             return Double.NEGATIVE_INFINITY;
  195.         }

  197.         final double xscale = x / scale;
  198.         final double logxscalepow = FastMath.log(xscale) * (shape - 1);

  200.         /*
  201.          * FastMath.pow(x / scale, shape) =
  202.          * FastMath.pow(xscale, shape) =
  203.          * FastMath.pow(xscale, shape - 1) * xscale
  204.          */
  205.         final double xscalepowshape = FastMath.exp(logxscalepow) * xscale;

  207.         return FastMath.log(shape / scale) + logxscalepow - xscalepowshape;
  208.     }

  210.     /** {@inheritDoc} */
  211.     public double cumulativeProbability(double x) {
  212.         double ret;
  213.         if (x <= 0.0) {
  214.             ret = 0.0;
  215.         } else {
  216.             ret = 1.0 - FastMath.exp(-FastMath.pow(x / scale, shape));
  217.         }
  218.         return ret;
  219.     }

  221.     /**
  222.      * {@inheritDoc}
  223.      *
  224.      * Returns {@code 0} when {@code p == 0} and
  225.      * {@code Double.POSITIVE_INFINITY} when {@code p == 1}.
  226.      */
  227.     @Override
  228.     public double inverseCumulativeProbability(double p) {
  229.         double ret;
  230.         if (p < 0.0 || p > 1.0) {
  231.             throw new OutOfRangeException(p, 0.0, 1.0);
  232.         } else if (p == 0) {
  233.             ret = 0.0;
  234.         } else  if (p == 1) {
  235.             ret = Double.POSITIVE_INFINITY;
  236.         } else {
  237.             ret = scale * FastMath.pow(-FastMath.log1p(-p), 1.0 / shape);
  238.         }
  239.         return ret;
  240.     }

  242.     /**
  243.      * Return the absolute accuracy setting of the solver used to estimate
  244.      * inverse cumulative probabilities.
  245.      *
  246.      * @return the solver absolute accuracy.
  247.      * @since 2.1
  248.      */
  249.     @Override
  250.     protected double getSolverAbsoluteAccuracy() {
  251.         return solverAbsoluteAccuracy;
  252.     }

  254.     /**
  255.      * {@inheritDoc}
  256.      *
  257.      * The mean is {@code scale * Gamma(1 + (1 / shape))}, where {@code Gamma()}
  258.      * is the Gamma-function.
  259.      */
  260.     public double getNumericalMean() {
  261.         if (!numericalMeanIsCalculated) {
  262.             numericalMean = calculateNumericalMean();
  263.             numericalMeanIsCalculated = true;
  264.         }
  265.         return numericalMean;
  266.     }

  268.     /**
  269.      * used by {@link #getNumericalMean()}
  270.      *
  271.      * @return the mean of this distribution
  272.      */
  273.     protected double calculateNumericalMean() {
  274.         final double sh = getShape();
  275.         final double sc = getScale();

  277.         return sc * FastMath.exp(Gamma.logGamma(1 + (1 / sh)));
  278.     }

  280.     /**
  281.      * {@inheritDoc}
  282.      *
  283.      * The variance is {@code scale^2 * Gamma(1 + (2 / shape)) - mean^2}
  284.      * where {@code Gamma()} is the Gamma-function.
  285.      */
  286.     public double getNumericalVariance() {
  287.         if (!numericalVarianceIsCalculated) {
  288.             numericalVariance = calculateNumericalVariance();
  289.             numericalVarianceIsCalculated = true;
  290.         }
  291.         return numericalVariance;
  292.     }

  294.     /**
  295.      * used by {@link #getNumericalVariance()}
  296.      *
  297.      * @return the variance of this distribution
  298.      */
  299.     protected double calculateNumericalVariance() {
  300.         final double sh = getShape();
  301.         final double sc = getScale();
  302.         final double mn = getNumericalMean();

  304.         return (sc * sc) * FastMath.exp(Gamma.logGamma(1 + (2 / sh))) -
  305.                (mn * mn);
  306.     }

  308.     /**
  309.      * {@inheritDoc}
  310.      *
  311.      * The lower bound of the support is always 0 no matter the parameters.
  312.      *
  313.      * @return lower bound of the support (always 0)
  314.      */
  315.     public double getSupportLowerBound() {
  316.         return 0;
  317.     }

  319.     /**
  320.      * {@inheritDoc}
  321.      *
  322.      * The upper bound of the support is always positive infinity
  323.      * no matter the parameters.
  324.      *
  325.      * @return upper bound of the support (always
  326.      * {@code Double.POSITIVE_INFINITY})
  327.      */
  328.     public double getSupportUpperBound() {
  329.         return Double.POSITIVE_INFINITY;
  330.     }

  332.     /** {@inheritDoc} */
  333.     public boolean isSupportLowerBoundInclusive() {
  334.         return true;
  335.     }

  337.     /** {@inheritDoc} */
  338.     public boolean isSupportUpperBoundInclusive() {
  339.         return false;
  340.     }

  342.     /**
  343.      * {@inheritDoc}
  344.      *
  345.      * The support of this distribution is connected.
  346.      *
  347.      * @return {@code true}
  348.      */
  349.     public boolean isSupportConnected() {
  350.         return true;
  351.     }
  352. }
Created with JaCoCo 0.7.2.201409121644