ROL
ROL_Gamma.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_GAMMA_HPP
11 #define ROL_GAMMA_HPP
12 
13 #include "ROL_Distribution.hpp"
14 #include "ROL_ParameterList.hpp"
15 
16 #include <math.h>
17 
18 namespace ROL {
19 
20 template<class Real>
21 class Gamma : public Distribution<Real> {
22 private:
23  Real shape_;
24  Real scale_;
26  Real coeff_;
27 
28  Real igamma(const Real s, const Real x) const {
29  Real sum = 0., term = 1./s;
30  size_t n = 1;
31  while ( term != 0. ) {
32  sum += term;
33  term *= x/(s+(Real)n);
34  n++;
35  }
36  return std::pow(x,s)*std::exp(-x)*sum;
37  }
38 
39 public:
40  Gamma(const Real shape = 1., const Real scale = 1.)
41  : shape_((shape > 0.) ? shape : 1.), scale_((scale > 0.) ? scale : 1.) {
42  gamma_shape_ = tgamma(shape_);
43  coeff_ = 1./(gamma_shape_*std::pow(scale_,shape_));
44  }
45 
46  Gamma(ROL::ParameterList &parlist) {
47  shape_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gamma").get("Shape",1.);
48  scale_ = parlist.sublist("SOL").sublist("Distribution").sublist("Gamma").get("Scale",1.);
49  shape_ = (shape_ > 0.) ? shape_ : 1.;
50  scale_ = (scale_ > 0.) ? scale_ : 1.;
51  gamma_shape_ = tgamma(shape_);
52  coeff_ = 1./(gamma_shape_*std::pow(scale_,shape_));
53  }
54 
55  Real evaluatePDF(const Real input) const {
56  return ((input <= 0.) ? 0. : coeff_*std::pow(input,shape_-1.)*std::exp(-input/scale_));
57  }
58 
59  Real evaluateCDF(const Real input) const {
60  return ((input <= 0.) ? 0. : igamma(shape_,input/scale_)/gamma_shape_);
61  }
62 
63  Real integrateCDF(const Real input) const {
64  Real x = input/scale_;
65  return ((input <= 0.) ? 0. : (x*igamma(shape_,x) - igamma(shape_+1.,x))/scale_);
66  }
67 
68  Real invertCDF(const Real input) const {
69  if ( input <= 0. ) {
70  return 0.;
71  }
72  Real x = input*gamma_shape_;
73  Real fx = evaluateCDF(x)-input;
74  Real s = 0., xs = 0., a = 1., tmp = 0.;
75  for (size_t i = 0; i < 100; i++) {
76  if ( std::abs(fx) < ROL_EPSILON<Real>() ) { break; }
77  s = -fx/evaluatePDF(x);
78  a = 1.0;
79  xs = x + a*s;
80  tmp = fx;
81  fx = evaluateCDF(xs)-input;
82  while ( std::abs(fx) > (1.0 - 1.e-4*a)*std::abs(tmp) ) {
83  a *= 0.5;
84  xs = x + a*s;
85  fx = evaluateCDF(xs)-input;
86  }
87  x = xs;
88  }
89  return x;
90  }
91 
92  Real moment(const size_t m) const {
93  if ( m == 1 ) {
94  return shape_*scale_;
95  }
96  if ( m == 2 ) {
97  return shape_*scale_*scale_*(1. + shape_);
98  }
99  return std::pow(scale_,m)*tgamma(shape_+(Real)m)/gamma_shape_;
100  }
101 
102  Real lowerBound(void) const {
103  return 0.;
104  }
105 
106  Real upperBound(void) const {
107  return ROL_INF<Real>();
108  }
109 
110  void test(std::ostream &outStream = std::cout ) const {
111  size_t size = 3;
112  std::vector<Real> X(size,0.);
113  std::vector<int> T(size,0);
114  X[0] = -4.0*(Real)rand()/(Real)RAND_MAX;
115  T[0] = 0;
116  X[1] = 0.;
117  T[1] = 1;
118  X[2] = 4.0*(Real)rand()/(Real)RAND_MAX;
119  T[2] = 0;
120  Distribution<Real>::test(X,T,outStream);
121  }
122 };
123 
124 }
125 
126 #endif
Real gamma_shape_
Definition: ROL_Gamma.hpp:25
Real invertCDF(const Real input) const
Definition: ROL_Gamma.hpp:68
Gamma(const Real shape=1., const Real scale=1.)
Definition: ROL_Gamma.hpp:40
Real igamma(const Real s, const Real x) const
Definition: ROL_Gamma.hpp:28
Gamma(ROL::ParameterList &parlist)
Definition: ROL_Gamma.hpp:46
Real evaluatePDF(const Real input) const
Definition: ROL_Gamma.hpp:55
Real evaluateCDF(const Real input) const
Definition: ROL_Gamma.hpp:59
void test(std::ostream &outStream=std::cout) const
Definition: ROL_Gamma.hpp:110
Real upperBound(void) const
Definition: ROL_Gamma.hpp:106
Real lowerBound(void) const
Definition: ROL_Gamma.hpp:102
Real scale_
Definition: ROL_Gamma.hpp:24
Real coeff_
Definition: ROL_Gamma.hpp:26
virtual void test(std::ostream &outStream=std::cout) const
Real shape_
Definition: ROL_Gamma.hpp:23
Real moment(const size_t m) const
Definition: ROL_Gamma.hpp:92
Real integrateCDF(const Real input) const
Definition: ROL_Gamma.hpp:63