ROL
ROL_Beta.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_BETA_HPP
45#define ROL_BETA_HPP
46
47#include "ROL_Distribution.hpp"
48#include "ROL_ParameterList.hpp"
49
50#include <math.h>
51
52namespace ROL {
53
54template<class Real>
55class Beta : public Distribution<Real> {
56private:
57 Real shape1_;
58 Real shape2_;
59 Real coeff_;
60
61 std::vector<Real> pts_;
62 std::vector<Real> wts_;
63
65 pts_.clear(); pts_.resize(20,0.); wts_.clear(); wts_.resize(20,0.);
66 wts_[0] = 0.1527533871307258; pts_[0] = -0.0765265211334973;
67 wts_[1] = 0.1527533871307258; pts_[1] = 0.0765265211334973;
68 wts_[2] = 0.1491729864726037; pts_[2] = -0.2277858511416451;
69 wts_[3] = 0.1491729864726037; pts_[3] = 0.2277858511416451;
70 wts_[4] = 0.1420961093183820; pts_[4] = -0.3737060887154195;
71 wts_[5] = 0.1420961093183820; pts_[5] = 0.3737060887154195;
72 wts_[6] = 0.1316886384491766; pts_[6] = -0.5108670019508271;
73 wts_[7] = 0.1316886384491766; pts_[7] = 0.5108670019508271;
74 wts_[8] = 0.1181945319615184; pts_[8] = -0.6360536807265150;
75 wts_[9] = 0.1181945319615184; pts_[9] = 0.6360536807265150;
76 wts_[10] = 0.1019301198172404; pts_[10] = -0.7463319064601508;
77 wts_[11] = 0.1019301198172404; pts_[11] = 0.7463319064601508;
78 wts_[12] = 0.0832767415767048; pts_[12] = -0.8391169718222188;
79 wts_[13] = 0.0832767415767048; pts_[13] = 0.8391169718222188;
80 wts_[14] = 0.0626720483341091; pts_[14] = -0.9122344282513259;
81 wts_[15] = 0.0626720483341091; pts_[15] = 0.9122344282513259;
82 wts_[16] = 0.0406014298003869; pts_[16] = -0.9639719272779138;
83 wts_[17] = 0.0406014298003869; pts_[17] = 0.9639719272779138;
84 wts_[18] = 0.0176140071391521; pts_[18] = -0.9931285991850949;
85 wts_[19] = 0.0176140071391521; pts_[19] = 0.9931285991850949;
86 for (size_t i = 0; i < 20; i++) {
87 wts_[i] *= 0.5;
88 pts_[i] += 1.; pts_[i] *= 0.5;
89 }
90 }
91
92 Real ibeta(const Real x) const {
93 Real pt = 0., wt = 0., sum = 0.;
94 for (size_t i = 0; i < pts_.size(); i++) {
95 wt = x*wts_[i];
96 pt = x*pts_[i];
97 sum += wt*std::pow(pt,shape1_-1)*std::pow(1.-pt,shape2_-1);
98 }
99 return sum;
100 }
101
102public:
103 Beta(const Real shape1 = 2., const Real shape2 = 2.)
104 : shape1_((shape1 > 0.) ? shape1 : 2.), shape2_((shape2 > 0.) ? shape2 : 2.) {
105 coeff_ = tgamma(shape1_+shape2_)/(tgamma(shape1_)*tgamma(shape2_));
107 }
108
109 Beta(ROL::ParameterList &parlist) {
110 shape1_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 1",2.);
111 shape2_ = parlist.sublist("SOL").sublist("Distribution").sublist("Beta").get("Shape 2",2.);
112 shape1_ = (shape1_ > 0.) ? shape1_ : 2.;
113 shape2_ = (shape2_ > 0.) ? shape2_ : 2.;
114 coeff_ = tgamma(shape1_+shape2_)/(tgamma(shape1_)*tgamma(shape2_));
116 }
117
118 Real evaluatePDF(const Real input) const {
119 return ((input > 0.) ? ((input < 1.) ?
120 coeff_*std::pow(input,shape1_-1.)*std::pow(1.-input,shape2_-1) : 0.) : 0.);
121 }
122
123 Real evaluateCDF(const Real input) const {
124 return ((input > 0.) ? ((input < 1.) ? coeff_*ibeta(input) : 1.) : 0.);
125 }
126
127 Real integrateCDF(const Real input) const {
128 ROL_TEST_FOR_EXCEPTION( true, std::invalid_argument,
129 ">>> ERROR (ROL::Beta): Beta integrateCDF not implemented!");
130 }
131
132 Real invertCDF(const Real input) const {
133 if ( input <= ROL_EPSILON<Real>() ) {
134 return 0.;
135 }
136 if ( input >= 1.-ROL_EPSILON<Real>() ) {
137 return 1.;
138 }
139 Real a = ROL_EPSILON<Real>(), b = 1.-ROL_EPSILON<Real>(), c = 0.;
140 Real fa = evaluateCDF(a) - input;
141 Real fc = 0.;
142 Real sa = ((fa < 0.) ? -1. : ((fa > 0.) ? 1. : 0.));
143 Real sc = 0.;
144 for (size_t i = 0; i < 100; i++) {
145 c = (a+b)*0.5;
146 fc = evaluateCDF(c) - input;
147 sc = ((fc < 0.) ? -1. : ((fc > 0.) ? 1. : 0.));
148 if ( fc == 0. || (b-a)*0.5 < ROL_EPSILON<Real>() ) {
149 break;
150 }
151 if ( sc == sa ) { a = c; fa = fc; sa = sc; }
152 else { b = c; }
153 }
154 return c;
155 }
156
157 Real moment(const size_t m) const {
158 if ( m == 1 ) {
159 return shape1_/(shape1_ + shape2_);
160 }
161 if ( m == 2 ) {
162 return shape1_*(shape2_/(shape1_+shape2_+1.) + shape1_)/std::pow(shape1_+shape2_,2);
163 }
164 Real val = 1.;
165 for (size_t i = 0; i < m; i++) {
166 val *= (shape1_ + (Real)i)/(shape1_ + shape2_ + (Real)i);
167 }
168 return val;
169 }
170
171 Real lowerBound(void) const {
172 return 0.;
173 }
174
175 Real upperBound(void) const {
176 return 1.;
177 }
178
179 void test(std::ostream &outStream = std::cout ) const {
180 size_t size = 5;
181 std::vector<Real> X(size,0.);
182 std::vector<int> T(size,0);
183 X[0] = -4.*(Real)rand()/(Real)RAND_MAX;
184 T[0] = 0;
185 X[1] = 0.;
186 T[1] = 1;
187 X[2] = 0.5*(Real)rand()/(Real)RAND_MAX;
188 T[2] = 0;
189 X[3] = 1.;
190 T[3] = 1;
191 X[4] = 1.+4.*(Real)rand()/(Real)RAND_MAX;
192 T[4] = 0;
193 Distribution<Real>::test(X,T,outStream);
194 }
195};
196
197}
198
199#endif
Real upperBound(void) const
Definition: ROL_Beta.hpp:175
Real shape1_
Definition: ROL_Beta.hpp:57
void initializeQuadrature(void)
Definition: ROL_Beta.hpp:64
std::vector< Real > pts_
Definition: ROL_Beta.hpp:61
Real lowerBound(void) const
Definition: ROL_Beta.hpp:171
Real moment(const size_t m) const
Definition: ROL_Beta.hpp:157
Real evaluateCDF(const Real input) const
Definition: ROL_Beta.hpp:123
Real evaluatePDF(const Real input) const
Definition: ROL_Beta.hpp:118
Real coeff_
Definition: ROL_Beta.hpp:59
Real shape2_
Definition: ROL_Beta.hpp:58
Real invertCDF(const Real input) const
Definition: ROL_Beta.hpp:132
Real integrateCDF(const Real input) const
Definition: ROL_Beta.hpp:127
void test(std::ostream &outStream=std::cout) const
Definition: ROL_Beta.hpp:179
Real ibeta(const Real x) const
Definition: ROL_Beta.hpp:92
Beta(ROL::ParameterList &parlist)
Definition: ROL_Beta.hpp:109
std::vector< Real > wts_
Definition: ROL_Beta.hpp:62
Beta(const Real shape1=2., const Real shape2=2.)
Definition: ROL_Beta.hpp:103
virtual void test(std::ostream &outStream=std::cout) const