488 lines · plain
1// Copyright John Maddock 2006.2// Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_STATS_WEIBULL_HPP8#define BOOST_STATS_WEIBULL_HPP9 10// http://www.itl.nist.gov/div898/handbook/eda/section3/eda3668.htm11// http://mathworld.wolfram.com/WeibullDistribution.html12 13#include <boost/math/tools/config.hpp>14#include <boost/math/tools/numeric_limits.hpp>15#include <boost/math/tools/type_traits.hpp>16#include <boost/math/tools/cstdint.hpp>17#include <boost/math/distributions/fwd.hpp>18#include <boost/math/special_functions/gamma.hpp>19#include <boost/math/special_functions/log1p.hpp>20#include <boost/math/special_functions/expm1.hpp>21#include <boost/math/distributions/detail/common_error_handling.hpp>22#include <boost/math/distributions/complement.hpp>23 24namespace boost{ namespace math25{26namespace detail{27 28template <class RealType, class Policy>29BOOST_MATH_GPU_ENABLED inline bool check_weibull_shape(30 const char* function,31 RealType shape,32 RealType* result, const Policy& pol)33{34 if((shape <= 0) || !(boost::math::isfinite)(shape))35 {36 *result = policies::raise_domain_error<RealType>(37 function,38 "Shape parameter is %1%, but must be > 0 !", shape, pol);39 return false;40 }41 return true;42}43 44template <class RealType, class Policy>45BOOST_MATH_GPU_ENABLED inline bool check_weibull_x(46 const char* function,47 RealType const& x,48 RealType* result, const Policy& pol)49{50 if((x < 0) || !(boost::math::isfinite)(x))51 {52 *result = policies::raise_domain_error<RealType>(53 function,54 "Random variate is %1% but must be >= 0 !", x, pol);55 return false;56 }57 return true;58}59 60template <class RealType, class Policy>61BOOST_MATH_GPU_ENABLED inline bool check_weibull(62 const char* function,63 RealType scale,64 RealType shape,65 RealType* result, const Policy& pol)66{67 return check_scale(function, scale, result, pol) && check_weibull_shape(function, shape, result, pol);68}69 70} // namespace detail71 72template <class RealType = double, class Policy = policies::policy<> >73class weibull_distribution74{75public:76 using value_type = RealType;77 using policy_type = Policy;78 79 BOOST_MATH_GPU_ENABLED explicit weibull_distribution(RealType l_shape, RealType l_scale = 1)80 : m_shape(l_shape), m_scale(l_scale)81 {82 RealType result;83 detail::check_weibull("boost::math::weibull_distribution<%1%>::weibull_distribution", l_scale, l_shape, &result, Policy());84 }85 86 BOOST_MATH_GPU_ENABLED RealType shape()const87 {88 return m_shape;89 }90 91 BOOST_MATH_GPU_ENABLED RealType scale()const92 {93 return m_scale;94 }95private:96 //97 // Data members:98 //99 RealType m_shape; // distribution shape100 RealType m_scale; // distribution scale101};102 103using weibull = weibull_distribution<double>;104 105#ifdef __cpp_deduction_guides106template <class RealType>107weibull_distribution(RealType)->weibull_distribution<typename boost::math::tools::promote_args<RealType>::type>;108template <class RealType>109weibull_distribution(RealType,RealType)->weibull_distribution<typename boost::math::tools::promote_args<RealType>::type>;110#endif111 112template <class RealType, class Policy>113BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> range(const weibull_distribution<RealType, Policy>& /*dist*/)114{ // Range of permissible values for random variable x.115 using boost::math::tools::max_value;116 return boost::math::pair<RealType, RealType>(static_cast<RealType>(0), max_value<RealType>());117}118 119template <class RealType, class Policy>120BOOST_MATH_GPU_ENABLED inline boost::math::pair<RealType, RealType> support(const weibull_distribution<RealType, Policy>& /*dist*/)121{ // Range of supported values for random variable x.122 // This is range where cdf rises from 0 to 1, and outside it, the pdf is zero.123 using boost::math::tools::max_value;124 using boost::math::tools::min_value;125 return boost::math::pair<RealType, RealType>(min_value<RealType>(), max_value<RealType>());126 // A discontinuity at x == 0, so only support down to min_value.127}128 129template <class RealType, class Policy>130BOOST_MATH_GPU_ENABLED inline RealType pdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x)131{132 BOOST_MATH_STD_USING // for ADL of std functions133 134 constexpr auto function = "boost::math::pdf(const weibull_distribution<%1%>, %1%)";135 136 RealType shape = dist.shape();137 RealType scale = dist.scale();138 139 RealType result = 0;140 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))141 return result;142 if(false == detail::check_weibull_x(function, x, &result, Policy()))143 return result;144 145 if(x == 0)146 {147 if(shape == 1)148 {149 return 1 / scale;150 }151 if(shape > 1)152 {153 return 0;154 }155 return policies::raise_overflow_error<RealType>(function, 0, Policy());156 }157 result = exp(-pow(x / scale, shape));158 result *= pow(x / scale, shape - 1) * shape / scale;159 160 return result;161}162 163template <class RealType, class Policy>164BOOST_MATH_GPU_ENABLED inline RealType logpdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x)165{166 BOOST_MATH_STD_USING // for ADL of std functions167 168 constexpr auto function = "boost::math::logpdf(const weibull_distribution<%1%>, %1%)";169 170 RealType shape = dist.shape();171 RealType scale = dist.scale();172 173 RealType result = 0;174 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))175 return result;176 if(false == detail::check_weibull_x(function, x, &result, Policy()))177 return result;178 179 if(x == 0)180 {181 if(shape == 1)182 {183 return log(1 / scale);184 }185 if(shape > 1)186 {187 return 0;188 }189 return policies::raise_overflow_error<RealType>(function, 0, Policy());190 }191 192 result = log(shape) - shape * log(scale) + log(x) * (shape - 1) - pow(x / scale, shape);193 194 return result;195}196 197template <class RealType, class Policy>198BOOST_MATH_GPU_ENABLED inline RealType cdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x)199{200 BOOST_MATH_STD_USING // for ADL of std functions201 202 constexpr auto function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)";203 204 RealType shape = dist.shape();205 RealType scale = dist.scale();206 207 RealType result = 0;208 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))209 return result;210 if(false == detail::check_weibull_x(function, x, &result, Policy()))211 return result;212 213 result = -boost::math::expm1(-pow(x / scale, shape), Policy());214 215 return result;216}217 218template <class RealType, class Policy>219BOOST_MATH_GPU_ENABLED inline RealType logcdf(const weibull_distribution<RealType, Policy>& dist, const RealType& x)220{221 BOOST_MATH_STD_USING // for ADL of std functions222 223 constexpr auto function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)";224 225 RealType shape = dist.shape();226 RealType scale = dist.scale();227 228 RealType result = 0;229 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))230 return result;231 if(false == detail::check_weibull_x(function, x, &result, Policy()))232 return result;233 234 result = log1p(-exp(-pow(x / scale, shape)), Policy());235 236 return result;237}238 239template <class RealType, class Policy>240BOOST_MATH_GPU_ENABLED inline RealType quantile(const weibull_distribution<RealType, Policy>& dist, const RealType& p)241{242 BOOST_MATH_STD_USING // for ADL of std functions243 244 constexpr auto function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)";245 246 RealType shape = dist.shape();247 RealType scale = dist.scale();248 249 RealType result = 0;250 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))251 return result;252 if(false == detail::check_probability(function, p, &result, Policy()))253 return result;254 255 if(p == 1)256 return policies::raise_overflow_error<RealType>(function, 0, Policy());257 258 result = scale * pow(-boost::math::log1p(-p, Policy()), 1 / shape);259 260 return result;261}262 263template <class RealType, class Policy>264BOOST_MATH_GPU_ENABLED inline RealType cdf(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c)265{266 BOOST_MATH_STD_USING // for ADL of std functions267 268 constexpr auto function = "boost::math::cdf(const weibull_distribution<%1%>, %1%)";269 270 RealType shape = c.dist.shape();271 RealType scale = c.dist.scale();272 273 RealType result = 0;274 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))275 return result;276 if(false == detail::check_weibull_x(function, c.param, &result, Policy()))277 return result;278 279 result = exp(-pow(c.param / scale, shape));280 281 return result;282}283 284template <class RealType, class Policy>285BOOST_MATH_GPU_ENABLED inline RealType logcdf(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c)286{287 BOOST_MATH_STD_USING // for ADL of std functions288 289 constexpr auto function = "boost::math::logcdf(const weibull_distribution<%1%>, %1%)";290 291 RealType shape = c.dist.shape();292 RealType scale = c.dist.scale();293 294 RealType result = 0;295 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))296 return result;297 if(false == detail::check_weibull_x(function, c.param, &result, Policy()))298 return result;299 300 result = -pow(c.param / scale, shape);301 302 return result;303}304 305template <class RealType, class Policy>306BOOST_MATH_GPU_ENABLED inline RealType quantile(const complemented2_type<weibull_distribution<RealType, Policy>, RealType>& c)307{308 BOOST_MATH_STD_USING // for ADL of std functions309 310 constexpr auto function = "boost::math::quantile(const weibull_distribution<%1%>, %1%)";311 312 RealType shape = c.dist.shape();313 RealType scale = c.dist.scale();314 RealType q = c.param;315 316 RealType result = 0;317 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))318 return result;319 if(false == detail::check_probability(function, q, &result, Policy()))320 return result;321 322 if(q == 0)323 return policies::raise_overflow_error<RealType>(function, 0, Policy());324 325 result = scale * pow(-log(q), 1 / shape);326 327 return result;328}329 330template <class RealType, class Policy>331BOOST_MATH_GPU_ENABLED inline RealType mean(const weibull_distribution<RealType, Policy>& dist)332{333 BOOST_MATH_STD_USING // for ADL of std functions334 335 constexpr auto function = "boost::math::mean(const weibull_distribution<%1%>)";336 337 RealType shape = dist.shape();338 RealType scale = dist.scale();339 340 RealType result = 0;341 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))342 return result;343 344 result = scale * boost::math::tgamma(1 + 1 / shape, Policy());345 return result;346}347 348template <class RealType, class Policy>349BOOST_MATH_GPU_ENABLED inline RealType variance(const weibull_distribution<RealType, Policy>& dist)350{351 RealType shape = dist.shape();352 RealType scale = dist.scale();353 354 constexpr auto function = "boost::math::variance(const weibull_distribution<%1%>)";355 356 RealType result = 0;357 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))358 {359 return result;360 }361 result = boost::math::tgamma(1 + 1 / shape, Policy());362 result *= -result;363 result += boost::math::tgamma(1 + 2 / shape, Policy());364 result *= scale * scale;365 return result;366}367 368template <class RealType, class Policy>369BOOST_MATH_GPU_ENABLED inline RealType mode(const weibull_distribution<RealType, Policy>& dist)370{371 BOOST_MATH_STD_USING // for ADL of std function pow.372 373 constexpr auto function = "boost::math::mode(const weibull_distribution<%1%>)";374 375 RealType shape = dist.shape();376 RealType scale = dist.scale();377 378 RealType result = 0;379 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))380 {381 return result;382 }383 if(shape <= 1)384 return 0;385 result = scale * pow((shape - 1) / shape, 1 / shape);386 return result;387}388 389template <class RealType, class Policy>390BOOST_MATH_GPU_ENABLED inline RealType median(const weibull_distribution<RealType, Policy>& dist)391{392 BOOST_MATH_STD_USING // for ADL of std function pow.393 394 constexpr auto function = "boost::math::median(const weibull_distribution<%1%>)";395 396 RealType shape = dist.shape(); // Wikipedia k397 RealType scale = dist.scale(); // Wikipedia lambda398 399 RealType result = 0;400 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))401 {402 return result;403 }404 using boost::math::constants::ln_two;405 result = scale * pow(ln_two<RealType>(), 1 / shape);406 return result;407}408 409template <class RealType, class Policy>410BOOST_MATH_GPU_ENABLED inline RealType skewness(const weibull_distribution<RealType, Policy>& dist)411{412 BOOST_MATH_STD_USING // for ADL of std functions413 414 constexpr auto function = "boost::math::skewness(const weibull_distribution<%1%>)";415 416 RealType shape = dist.shape();417 RealType scale = dist.scale();418 419 RealType result = 0;420 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))421 {422 return result;423 }424 425 RealType g1 = boost::math::tgamma(1 + 1 / shape, Policy());426 RealType g2 = boost::math::tgamma(1 + 2 / shape, Policy());427 RealType g3 = boost::math::tgamma(1 + 3 / shape, Policy());428 RealType d = pow(g2 - g1 * g1, RealType(1.5));429 430 result = (2 * g1 * g1 * g1 - 3 * g1 * g2 + g3) / d;431 return result;432}433 434template <class RealType, class Policy>435BOOST_MATH_GPU_ENABLED inline RealType kurtosis_excess(const weibull_distribution<RealType, Policy>& dist)436{437 BOOST_MATH_STD_USING // for ADL of std functions438 439 constexpr auto function = "boost::math::kurtosis_excess(const weibull_distribution<%1%>)";440 441 RealType shape = dist.shape();442 RealType scale = dist.scale();443 444 RealType result = 0;445 if(false == detail::check_weibull(function, scale, shape, &result, Policy()))446 return result;447 448 RealType g1 = boost::math::tgamma(1 + 1 / shape, Policy());449 RealType g2 = boost::math::tgamma(1 + 2 / shape, Policy());450 RealType g3 = boost::math::tgamma(1 + 3 / shape, Policy());451 RealType g4 = boost::math::tgamma(1 + 4 / shape, Policy());452 RealType g1_2 = g1 * g1;453 RealType g1_4 = g1_2 * g1_2;454 RealType d = g2 - g1_2;455 d *= d;456 457 result = -6 * g1_4 + 12 * g1_2 * g2 - 3 * g2 * g2 - 4 * g1 * g3 + g4;458 result /= d;459 return result;460}461 462template <class RealType, class Policy>463BOOST_MATH_GPU_ENABLED inline RealType kurtosis(const weibull_distribution<RealType, Policy>& dist)464{465 return kurtosis_excess(dist) + 3;466}467 468template <class RealType, class Policy>469BOOST_MATH_GPU_ENABLED inline RealType entropy(const weibull_distribution<RealType, Policy>& dist)470{471 BOOST_MATH_STD_USING472 RealType k = dist.shape();473 RealType lambda = dist.scale();474 return constants::euler<RealType>()*(1-1/k) + log(lambda/k) + 1;475}476 477} // namespace math478} // namespace boost479 480// This include must be at the end, *after* the accessors481// for this distribution have been defined, in order to482// keep compilers that support two-phase lookup happy.483#include <boost/math/distributions/detail/derived_accessors.hpp>484 485#endif // BOOST_STATS_WEIBULL_HPP486 487 488