Line data    Source code

       1             : //* This file is part of the MOOSE framework
       2             : //* https://mooseframework.inl.gov
       3             : //*
       4             : //* All rights reserved, see COPYRIGHT for full restrictions
       5             : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
       6             : //*
       7             : //* Licensed under LGPL 2.1, please see LICENSE for details
       8             : //* https://www.gnu.org/licenses/lgpl-2.1.html
       9             : 
      10             : #include "Gamma.h"
      11             : #include "math.h"
      12             : #include "libmesh/utility.h"
      13             : 
      14             : registerMooseObject("StochasticToolsApp", Gamma);
      15             : 
      16             : InputParameters
      17          64 : Gamma::validParams()
      18             : {
      19          64 :   InputParameters params = Distribution::validParams();
      20          64 :   params.addClassDescription("Gamma distribution");
      21         128 :   params.addRequiredRangeCheckedParam<Real>("shape", "shape > 0", "Shape parameter (k or alpha).");
      22         192 :   params.addRangeCheckedParam<Real>(
      23         128 :       "scale", 1.0, "scale > 0", "Scale parameter (theta or 1/beta).");
      24          64 :   return params;
      25           0 : }
      26             : 
      27          32 : Gamma::Gamma(const InputParameters & parameters)
      28          96 :   : Distribution(parameters), _alpha(getParam<Real>("shape")), _theta(getParam<Real>("scale"))
      29             : {
      30          32 : }
      31             : 
      32             : Real
      33          32 : Gamma::pdf(const Real & x, const Real & alpha, const Real & beta)
      34             : {
      35          32 :   if (x <= 0.0)
      36             :     return 0.0;
      37             :   else
      38          32 :     return std::pow(beta, alpha) * std::pow(x, alpha - 1.0) * std::exp(-beta * x) /
      39          32 :            std::tgamma(alpha);
      40             : }
      41             : 
      42             : Real
      43          32 : Gamma::cdf(const Real & x, const Real & alpha, const Real & beta)
      44             : {
      45          32 :   if (x <= 0.0)
      46             :     return 0.0;
      47             :   else
      48          32 :     return incompleteGamma(alpha, beta * x);
      49             : }
      50             : 
      51             : Real
      52          32 : Gamma::quantile(const Real & p, const Real & alpha, const Real & beta)
      53             : {
      54          32 :   return incompleteGammaInv(alpha, p) / beta;
      55             : }
      56             : 
      57             : Real
      58          32 : Gamma::pdf(const Real & x) const
      59             : {
      60          32 :   return pdf(x, _alpha, 1.0 / _theta);
      61             : }
      62             : 
      63             : Real
      64          32 : Gamma::cdf(const Real & x) const
      65             : {
      66          32 :   return cdf(x, _alpha, 1.0 / _theta);
      67             : }
      68             : 
      69             : Real
      70          32 : Gamma::quantile(const Real & p) const
      71             : {
      72          32 :   return quantile(p, _alpha, 1.0 / _theta);
      73             : }
      74             : 
      75             : Real
      76         192 : Gamma::incompleteGamma(const Real & a, const Real & x)
      77             : {
      78             :   const Real tol = 1e-14;
      79             :   const unsigned int max_iter = 1e6;
      80         192 :   const Real coef = std::pow(x, a) * std::exp(-x) / std::tgamma(a + 1.0);
      81             :   Real cn = 1.0;
      82             :   Real an = a;
      83             :   Real val = 1.0;
      84        3984 :   for (unsigned int i = 0; i < max_iter; ++i)
      85             :   {
      86        3984 :     an += 1.0;
      87        3984 :     cn *= x / an;
      88        3984 :     val += cn;
      89        3984 :     if (std::abs(cn / val) < tol)
      90         192 :       return coef * val;
      91             :   }
      92             : 
      93             :   mooseAssert(false, "Could not compute incomplete gamma function.");
      94           0 :   return coef * val;
      95             : }
      96             : 
      97             : Real
      98          32 : Gamma::incompleteGammaInv(const Real & a, const Real & p)
      99             : {
     100          32 :   Real x = a > 1.0 ? a : std::pow(p * std::tgamma(a + 1.0), 1.0 / a);
     101          32 :   const Real scale = std::tgamma(a);
     102             :   const Real tol = 1e-14;
     103             :   const unsigned int max_iter = 1e6;
     104             :   Real f;
     105             :   Real df;
     106         160 :   for (unsigned int i = 0; i < max_iter; ++i)
     107             :   {
     108         160 :     f = incompleteGamma(a, x) - p;
     109         160 :     if (std::abs(f) < tol)
     110          32 :       return x;
     111         128 :     df = std::pow(x, a - 1.0) * std::exp(-x) / scale;
     112         128 :     x -= f / df;
     113             :   }
     114             : 
     115             :   mooseAssert(false, "Could not find inverse of incomplete gamma function.");
     116           0 :   return x;
     117             : }