ROL
ROL_Beta.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_BETA_HPP
11 #define ROL_BETA_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 Beta : public Distribution<Real> {
22 private:
23  Real shape1_;
24  Real shape2_;
25  Real coeff_;
26 
27  std::vector<Real> pts_;
28  std::vector<Real> wts_;
29 
30  void initializeQuadrature(void) {
31  pts_.clear(); pts_.resize(20,0.); wts_.clear(); wts_.resize(20,0.);
32  wts_[0] = 0.1527533871307258; pts_[0] = -0.0765265211334973;
33  wts_[1] = 0.1527533871307258; pts_[1] = 0.0765265211334973;
34  wts_[2] = 0.1491729864726037; pts_[2] = -0.2277858511416451;
35  wts_[3] = 0.1491729864726037; pts_[3] = 0.2277858511416451;
36  wts_[4] = 0.1420961093183820; pts_[4] = -0.3737060887154195;
37  wts_[5] = 0.1420961093183820; pts_[5] = 0.3737060887154195;
38  wts_[6] = 0.1316886384491766; pts_[6] = -0.5108670019508271;
39  wts_[7] = 0.1316886384491766; pts_[7] = 0.5108670019508271;
40  wts_[8] = 0.1181945319615184; pts_[8] = -0.6360536807265150;
41  wts_[9] = 0.1181945319615184; pts_[9] = 0.6360536807265150;
42  wts_[10] = 0.1019301198172404; pts_[10] = -0.7463319064601508;
43  wts_[11] = 0.1019301198172404; pts_[11] = 0.7463319064601508;
44  wts_[12] = 0.0832767415767048; pts_[12] = -0.8391169718222188;
45  wts_[13] = 0.0832767415767048; pts_[13] = 0.8391169718222188;
46  wts_[14] = 0.0626720483341091; pts_[14] = -0.9122344282513259;
47  wts_[15] = 0.0626720483341091; pts_[15] = 0.9122344282513259;
48  wts_[16] = 0.0406014298003869; pts_[16] = -0.9639719272779138;
49  wts_[17] = 0.0406014298003869; pts_[17] = 0.9639719272779138;
50  wts_[18] = 0.0176140071391521; pts_[18] = -0.9931285991850949;
51  wts_[19] = 0.0176140071391521; pts_[19] = 0.9931285991850949;
52  for (size_t i = 0; i < 20; i++) {
53  wts_[i] *= 0.5;
54  pts_[i] += 1.; pts_[i] *= 0.5;
55  }
56  }
57 
58  Real ibeta(const Real x) const {
59  Real pt = 0., wt = 0., sum = 0.;
60  for (size_t i = 0; i < pts_.size(); i++) {
61  wt = x*wts_[i];
62  pt = x*pts_[i];
63  sum += wt*std::pow(pt,shape1_-1)*std::pow(1.-pt,shape2_-1);
64  }
65  return sum;
66  }
67 
68 public:
69  Beta(const Real shape1 = 2., const Real shape2 = 2.)
70  : shape1_((shape1 > 0.) ? shape1 : 2.), shape2_((shape2 > 0.) ? shape2 : 2.) {
71  coeff_ = tgamma(shape1_+shape2_)/(tgamma(shape1_)*tgamma(shape2_));
73  }
74 
75  Beta(ROL::ParameterList &parlist) {
76  shape1_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 1",2.);
77  shape2_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 2",2.);
78  shape1_ = (shape1_ > 0.) ? shape1_ : 2.;
79  shape2_ = (shape2_ > 0.) ? shape2_ : 2.;
80  coeff_ = tgamma(shape1_+shape2_)/(tgamma(shape1_)*tgamma(shape2_));
82  }
83 
84  Real evaluatePDF(const Real input) const {
85  return ((input > 0.) ? ((input < 1.) ?
86  coeff_*std::pow(input,shape1_-1.)*std::pow(1.-input,shape2_-1) : 0.) : 0.);
87  }
88 
89  Real evaluateCDF(const Real input) const {
90  return ((input > 0.) ? ((input < 1.) ? coeff_*ibeta(input) : 1.) : 0.);
91  }
92 
93  Real integrateCDF(const Real input) const {
94  ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
95  ">>> ERROR (ROL::Beta): Beta integrateCDF not implemented!");
96  }
97 
98  Real invertCDF(const Real input) const {
99  if ( input <= ROL_EPSILON<Real>() ) {
100  return 0.;
101  }
102  if ( input >= 1.-ROL_EPSILON<Real>() ) {
103  return 1.;
104  }
105  Real a = ROL_EPSILON<Real>(), b = 1.-ROL_EPSILON<Real>(), c = 0.;
106  Real fa = evaluateCDF(a) - input;
107  Real fc = 0.;
108  Real sa = ((fa < 0.) ? -1. : ((fa > 0.) ? 1. : 0.));
109  Real sc = 0.;
110  for (size_t i = 0; i < 100; i++) {
111  c = (a+b)*0.5;
112  fc = evaluateCDF(c) - input;
113  sc = ((fc < 0.) ? -1. : ((fc > 0.) ? 1. : 0.));
114  if ( fc == 0. || (b-a)*0.5 < ROL_EPSILON<Real>() ) {
115  break;
116  }
117  if ( sc == sa ) { a = c; fa = fc; sa = sc; }
118  else { b = c; }
119  }
120  return c;
121  }
122 
123  Real moment(const size_t m) const {
124  if ( m == 1 ) {
125  return shape1_/(shape1_ + shape2_);
126  }
127  if ( m == 2 ) {
128  return shape1_*(shape2_/(shape1_+shape2_+1.) + shape1_)/std::pow(shape1_+shape2_,2);
129  }
130  Real val = 1.;
131  for (size_t i = 0; i < m; i++) {
132  val *= (shape1_ + (Real)i)/(shape1_ + shape2_ + (Real)i);
133  }
134  return val;
135  }
136 
137  Real lowerBound(void) const {
138  return 0.;
139  }
140 
141  Real upperBound(void) const {
142  return 1.;
143  }
144 
145  void test(std::ostream &outStream = std::cout ) const {
146  size_t size = 5;
147  std::vector<Real> X(size,0.);
148  std::vector<int> T(size,0);
149  X[0] = -4.*(Real)rand()/(Real)RAND_MAX;
150  T[0] = 0;
151  X[1] = 0.;
152  T[1] = 1;
153  X[2] = 0.5*(Real)rand()/(Real)RAND_MAX;
154  T[2] = 0;
155  X[3] = 1.;
156  T[3] = 1;
157  X[4] = 1.+4.*(Real)rand()/(Real)RAND_MAX;
158  T[4] = 0;
159  Distribution<Real>::test(X,T,outStream);
160  }
161 };
162 
163 }
164 
165 #endif
Real evaluatePDF(const Real input) const
Definition: ROL_Beta.hpp:84
Beta(const Real shape1=2., const Real shape2=2.)
Definition: ROL_Beta.hpp:69
Real integrateCDF(const Real input) const
Definition: ROL_Beta.hpp:93
Real shape2_
Definition: ROL_Beta.hpp:24
Real shape1_
Definition: ROL_Beta.hpp:23
Real ibeta(const Real x) const
Definition: ROL_Beta.hpp:58
void initializeQuadrature(void)
Definition: ROL_Beta.hpp:30
Real coeff_
Definition: ROL_Beta.hpp:25
Real lowerBound(void) const
Definition: ROL_Beta.hpp:137
void test(std::ostream &outStream=std::cout) const
Definition: ROL_Beta.hpp:145
Real upperBound(void) const
Definition: ROL_Beta.hpp:141
virtual void test(std::ostream &outStream=std::cout) const
std::vector< Real > pts_
Definition: ROL_Beta.hpp:27
Real evaluateCDF(const Real input) const
Definition: ROL_Beta.hpp:89
std::vector< Real > wts_
Definition: ROL_Beta.hpp:28
Beta(ROL::ParameterList &parlist)
Definition: ROL_Beta.hpp:75
Real moment(const size_t m) const
Definition: ROL_Beta.hpp:123
Real invertCDF(const Real input) const
Definition: ROL_Beta.hpp:98