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