1 // Arguments: Doubles, Doubles, Doubles 2 #include <stan/math/prim.hpp> 3 4 using stan::math::square; 5 using stan::math::var; 6 using std::numeric_limits; 7 using std::vector; 8 9 class AgradCdfLogNormal : public AgradCdfLogTest { 10 public: valid_values(vector<vector<double>> & parameters,vector<double> & cdf_log)11 void valid_values(vector<vector<double> >& parameters, 12 vector<double>& cdf_log) { 13 vector<double> param(3); 14 15 param[0] = 0; // y 16 param[1] = 0; // mu 17 param[2] = 1; // sigma 18 parameters.push_back(param); 19 cdf_log.push_back(-0.6931471805599452862268); // expected cdf_log 20 21 param[0] = 1; // y 22 param[1] = 0; // mu 23 param[2] = 1; // sigma 24 parameters.push_back(param); 25 cdf_log.push_back(-0.1727537790234499326392); // expected cdf_log 26 27 param[0] = -2; // y 28 param[1] = 0; // mu 29 param[2] = 1.5; // sigma 30 parameters.push_back(param); 31 cdf_log.push_back(-2.3945773661586438052); // expected cdf_log 32 33 param[0] = -3.5; // y 34 param[1] = 1.9; // mu 35 param[2] = 7.2; // sigma 36 parameters.push_back(param); 37 cdf_log.push_back(-1.484448229919656192521); // expected cdf_log 38 39 param[0] = 7.5; // y 40 param[1] = 0; // mu 41 param[2] = 1; // sigma 42 parameters.push_back(param); 43 cdf_log.push_back(-3.1908916729109475e-14); // expected cdf_log 44 } 45 invalid_values(vector<size_t> & index,vector<double> & value)46 void invalid_values(vector<size_t>& index, vector<double>& value) { 47 // y 48 49 // mu 50 index.push_back(1U); 51 value.push_back(numeric_limits<double>::infinity()); 52 53 index.push_back(1U); 54 value.push_back(-numeric_limits<double>::infinity()); 55 56 // sigma 57 index.push_back(2U); 58 value.push_back(0.0); 59 60 index.push_back(2U); 61 value.push_back(-1.0); 62 63 index.push_back(2U); 64 value.push_back(-numeric_limits<double>::infinity()); 65 } 66 has_lower_bound()67 bool has_lower_bound() { return false; } 68 has_upper_bound()69 bool has_upper_bound() { return false; } 70 71 template <typename T_y, typename T_loc, typename T_scale, typename T3, 72 typename T4, typename T5> cdf_log(const T_y & y,const T_loc & mu,const T_scale & sigma,const T3 &,const T4 &,const T5 &)73 stan::return_type_t<T_y, T_loc, T_scale> cdf_log(const T_y& y, 74 const T_loc& mu, 75 const T_scale& sigma, 76 const T3&, const T4&, 77 const T5&) { 78 return stan::math::normal_cdf_log(y, mu, sigma); 79 } 80 81 template <typename T_y, typename T_loc, typename T_scale, typename T3, 82 typename T4, typename T5> cdf_log_function(const T_y & y,const T_loc & mu,const T_scale & sigma,const T3 &,const T4 &,const T5 &)83 stan::return_type_t<T_y, T_loc, T_scale> cdf_log_function( 84 const T_y& y, const T_loc& mu, const T_scale& sigma, const T3&, const T4&, 85 const T5&) { 86 using stan::math::LOG_HALF; 87 using stan::math::SQRT_PI; 88 using stan::math::SQRT_TWO; 89 using std::exp; 90 using std::log; 91 using std::log1p; 92 using std::pow; 93 94 stan::return_type_t<T_y, T_loc, T_scale> cdf_log(0.0); 95 96 // define constants used in numerical approximation of log(CDF) for x<-20 97 stan::return_type_t<T_y, T_loc, T_scale> p[6] 98 = {-0.000658749161529837803157, -0.0160837851487422766278, 99 -0.125781726111229246204, -0.360344899949804439429, 100 -0.305326634961232344035, -0.0163153871373020978498}; 101 stan::return_type_t<T_y, T_loc, T_scale> q[6] 102 = {0.00233520497626869185443, 0.0605183413124413191178, 103 0.527905102951428412248, 1.87295284992346047209, 104 2.56852019228982242072, 1.0}; 105 stan::return_type_t<T_y, T_loc, T_scale> scaled_diff(0.0); 106 scaled_diff += (y - mu) / (sigma * SQRT_TWO); 107 stan::return_type_t<T_y, T_loc, T_scale> sigma2pi(0.0); 108 sigma2pi += sigma * SQRT_TWO; 109 110 // use erfc() instead of erf() in order to retain precision since for x>0 111 // erfc()->0 112 if (scaled_diff > 0.0) { 113 // CDF(x) = 1/2 + 1/2erf(x) = 1 - 1/2erfc(x) 114 cdf_log += log1p(-0.5 * erfc(scaled_diff)); 115 } else if (scaled_diff > -20.0) { 116 // CDF(x) = 1/2 - 1/2erf(-x) = 1/2erfc(-x) 117 cdf_log += log(erfc(-scaled_diff)) + LOG_HALF; 118 } else { 119 // check whether scaled_diff^10 term will overflow 120 if (10.0 * log(-scaled_diff) < log(std::numeric_limits<double>::max())) { 121 // entering territory where erfc(-x)~0 122 // need to use direct numerical approximation of cdf_log instead 123 // the following based on W. J. Cody, Math. Comp. 23(107):631-638 (1969) 124 // CDF(x) = 1/2erfc(-x) 125 stan::return_type_t<T_y, T_loc, T_scale> temp_p = 0.0; 126 stan::return_type_t<T_y, T_loc, T_scale> temp_q = 0.0; 127 stan::return_type_t<T_y, T_loc, T_scale> x2 = square(scaled_diff); 128 stan::return_type_t<T_y, T_loc, T_scale> pow_x = 0.0; 129 for (int j = 0; j < 6; j++) { 130 pow_x = pow(scaled_diff, -2 * j); 131 temp_p = temp_p - p[j] * pow_x; 132 temp_q = temp_q - q[j] * pow_x; 133 } 134 cdf_log += log(0.5) + log(1.0 / SQRT_PI + (temp_p / temp_q) / x2) 135 - log(-scaled_diff) - x2; 136 } else { 137 cdf_log += stan::math::negative_infinity(); 138 } 139 } 140 return cdf_log; 141 } 142 }; 143