44 #ifndef ROL_GAUSSIAN_HPP 45 #define ROL_GAUSSIAN_HPP 48 #include "Teuchos_ParameterList.hpp" 63 Real
erfi(
const Real p)
const {
64 Real val = 0., z = 0.;
66 z = std::sqrt(-std::log((1.+p)*0.5));
67 val = -(((c_[3]*z+c_[2])*z+c_[1])*z + c_[0])/((d_[1]*z+d_[0])*z + 1.);
72 val = p*(((a_[3]*z+a_[2])*z+a_[1])*z + a_[0])/((((b_[3]*z+b_[2])*z+b_[1])*z+b_[0])*z+1.);
75 z = std::sqrt(-std::log((1.-p)*0.5));
76 val = (((c_[3]*z+c_[2])*z+c_[1])*z+c_[0])/((d_[1]*z+d_[0])*z+1.);
79 val -= (erf(val)-p)/(2.0/std::sqrt(M_PI) * std::exp(-val*val));
80 val -= (erf(val)-p)/(2.0/std::sqrt(M_PI) * std::exp(-val*val));
86 Gaussian(
const Real mean = 0.,
const Real variance = 1.)
87 : mean_(mean), variance_((variance>0.) ? variance : 1.) {
88 a_.clear(); a_.resize(4,0.); b_.clear(); b_.resize(4,0.);
89 c_.clear(); c_.resize(4,0.); d_.clear(); d_.resize(2,0.);
90 a_[0] = 0.886226899; a_[1] = -1.645349621; a_[2] = 0.914624893; a_[3] = -0.140543331;
91 b_[0] = -2.118377725; b_[1] = 1.442710462; b_[2] = -0.329097515; b_[3] = 0.012229801;
92 c_[0] = -1.970840454; c_[1] = -1.624906493; c_[2] = 3.429567803; c_[3] = 1.641345311;
93 d_[0] = 3.543889200; d_[1] = 1.637067800;
97 mean_ = parlist.sublist(
"SOL").sublist(
"Distribution").sublist(
"Gaussian").get(
"Mean",0.);
98 variance_ = parlist.sublist(
"SOL").sublist(
"Distribution").sublist(
"Gaussian").get(
"Variance",1.);
99 variance_ = (variance_ > 0.) ? variance_ : 1.;
100 a_.clear(); a_.resize(4,0.); b_.clear(); b_.resize(4,0.);
101 c_.clear(); c_.resize(4,0.); d_.clear(); d_.resize(2,0.);
102 a_[0] = 0.886226899; a_[1] = -1.645349621; a_[2] = 0.914624893; a_[3] = -0.140543331;
103 b_[0] = -2.118377725; b_[1] = 1.442710462; b_[2] = -0.329097515; b_[3] = 0.012229801;
104 c_[0] = -1.970840454; c_[1] = -1.624906493; c_[2] = 3.429567803; c_[3] = 1.641345311;
105 d_[0] = 3.543889200; d_[1] = 1.637067800;
109 return std::exp(-std::pow(input-mean_,2)/(2.*variance_))/(std::sqrt(2.*M_PI*variance_));
113 return 0.5*(1.+erf((input-mean_)/std::sqrt(2.*variance_)));
117 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
118 ">>> ERROR (ROL::Gaussian): Gaussian integrateCDF not implemented!");
119 return ((input < mean_) ? 0.0 : input);
123 return std::sqrt(2.*variance_)*
erfi(2.*input-1.) +
mean_;
129 case 1: val =
mean_;
break;
130 case 2: val = std::pow(mean_,2) +
variance_;
break;
131 case 3: val = std::pow(mean_,3)
133 case 4: val = std::pow(mean_,4)
134 + 6.*std::pow(mean_,2)*variance_
135 + 3.*std::pow(variance_,2);
break;
136 case 5: val = std::pow(mean_,5)
137 + 10.*std::pow(mean_,3)*variance_
138 + 15.*mean_*std::pow(variance_,2);
break;
139 case 6: val = std::pow(mean_,6)
140 + 15.*std::pow(mean_,4)*variance_
141 + 45.*std::pow(mean_*variance_,2)
142 + 15.*std::pow(variance_,3);
break;
143 case 7: val = std::pow(mean_,7)
144 + 21.*std::pow(mean_,5)*variance_
145 + 105.*std::pow(mean_,3)*std::pow(variance_,2)
146 + 105.*mean_*std::pow(variance_,3);
break;
147 case 8: val = std::pow(mean_,8)
148 + 28.*std::pow(mean_,6)*variance_
149 + 210.*std::pow(mean_,4)*std::pow(variance_,2)
150 + 420.*std::pow(mean_,2)*std::pow(variance_,3)
151 + 105.*std::pow(variance_,4);
break;
153 TEUCHOS_TEST_FOR_EXCEPTION(
true, std::invalid_argument,
154 ">>> ERROR (ROL::Distribution): Gaussian moment not implemented for m > 8!");
167 void test(std::ostream &outStream = std::cout )
const {
169 std::vector<Real> X(size,4.*(Real)rand()/(Real)RAND_MAX - 2.);
170 std::vector<int> T(size,0);
Gaussian(Teuchos::ParameterList &parlist)
Gaussian(const Real mean=0., const Real variance=1.)
Real lowerBound(void) const
Real integrateCDF(const Real input) const
Real evaluateCDF(const Real input) const
virtual void test(std::ostream &outStream=std::cout) const
static const double ROL_INF
Real moment(const size_t m) const
static const double ROL_NINF
Real invertCDF(const Real input) const
Real upperBound(void) const
Real erfi(const Real p) const
Real evaluatePDF(const Real input) const
void test(std::ostream &outStream=std::cout) const