1/***************************************************************************** 2 * * 3 * UNURAN -- Universal Non-Uniform Random number generator * 4 * * 5 ***************************************************************************** 6 * * 7 * FILE: ninv_newton.ch * 8 * * 9 * Routines for Newton's root finding algorithm. * 10 * * 11 ***************************************************************************** 12 * * 13 * Copyright (c) 2000-2009 Wolfgang Hoermann and Josef Leydold * 14 * Department of Statistics and Mathematics, WU Wien, Austria * 15 * * 16 ***************************************************************************** 17 * * 18 * REFERENCES: * 19 * [1] Neumaier A. (to be published): Introduction to numerical analysis, * 20 * Cambridge University Press * 21 * * 22 ***************************************************************************** 23 * * 24 * For a given U find X with F(X) - U = 0. * 25 * * 26 * Method: damped Newton method. * 27 * * 28 *****************************************************************************/ 29 30/*****************************************************************************/ 31/** Newton method **/ 32/*****************************************************************************/ 33 34/*---------------------------------------------------------------------------*/ 35/* Constants */ 36 37/* maximal number of steps to leave flat region */ 38#define MAX_FLAT_COUNT (40) 39 40/*---------------------------------------------------------------------------*/ 41 42double 43_unur_ninv_newton( const struct unur_gen *gen, double U ) 44 /*----------------------------------------------------------------------*/ 45 /* sample from generator (use Newton's method) */ 46 /* */ 47 /* parameters: */ 48 /* gen ... pointer to generator object */ 49 /* U ... random number (uniform distribution) */ 50 /* */ 51 /* return: */ 52 /* double (sample from random variate) */ 53 /* */ 54 /* error: */ 55 /* return INFINITY */ 56 /*----------------------------------------------------------------------*/ 57{ 58 double x; /* point for netwon-iteration */ 59 double fx; /* cdf at x */ 60 double dfx; /* pdf at x */ 61 double fxabs; /* absolute valuo of fx */ 62 double xtmp, fxtmp; /* temprary variables for x and fx */ 63 double xold; /* remember last values for stopping criterion */ 64 double fxtmpabs; /* fabs of fxtmp */ 65 double damp; /* damping factor */ 66 double step; /* helps to escape from flat regions of the cdf */ 67 int i; /* counter for for-loop, index */ 68 int flat_count; /* counter of steps in flat region */ 69 double rel_u_resolution; /* relative u resolution */ 70 int x_goal, u_goal; /* whether precision goal is reached */ 71 72 /* check arguments */ 73 CHECK_NULL(gen,INFINITY); COOKIE_CHECK(gen,CK_NINV_GEN,INFINITY); 74 75 /* compute relative u resolution */ 76 rel_u_resolution = ( (GEN->u_resolution > 0.) ? 77 (GEN->Umax - GEN->Umin) * GEN->u_resolution : 78 INFINITY ); 79 80 /* -- 1. initialize starting interval -- */ 81 82 if (GEN->table_on) { 83 /* -- 1a. use table -- */ 84 85 /* 0 <= i <= table_size-2 */ 86 if ( _unur_FP_same(GEN->CDFmin,GEN->CDFmax) ) { 87 /* CDF values in table too close, so we use median point since */ 88 /* there is no difference between CDF values. */ 89 i = GEN->table_size/2; 90 } 91 else { 92 i = (int) ( GEN->table_size * (U - GEN->CDFmin) / (GEN->CDFmax - GEN->CDFmin) ); 93 if (i<0) i = 0; 94 else if (i > GEN->table_size - 2) i = GEN->table_size - 2; 95 } 96 97 if (_unur_FP_is_infinity(GEN->table[i+1])) { 98 x = GEN->table[i]; 99 fx = GEN->f_table[i]; 100 } 101 else { 102 x = GEN->table[i+1]; 103 fx = GEN->f_table[i+1]; 104 } 105 106 } 107 108 else { 109 /* 1b. -- no table available -- */ 110 x = GEN->s[0]; 111 fx = GEN->CDFs[0]; 112 } 113 114 /* -- 1c. check for boundary of truncated domain -- */ 115 116 if ( x < DISTR.trunc[0] ){ 117 x = DISTR.trunc[0]; 118 fx = GEN->Umin; /* = CDF(x) */ 119 } 120 else if ( x > DISTR.trunc[1] ){ 121 x = DISTR.trunc[1]; 122 fx = GEN->Umax; /* = CDF(x) */ 123 } 124 125 /* -- 1z. compute values for starting point -- */ 126 127 fx -= U; 128 dfx = PDF(x); 129 fxabs = fabs(fx); 130 xold = x; /* there is no old value yet */ 131 132 damp = 2.; /* to be halved at least once */ 133 step = 1.; 134 135 /* -- 2. Newton iteration -- */ 136 137 for (i=0; i < GEN->max_iter; i++) { 138 139 flat_count = 0; 140 while (_unur_iszero(dfx)) { /* function flat at x */ 141 /* printf("step: %g, x: %g, fx: %g, dfx: %g\n",step, x, fx, dfx); */ 142 143 if (_unur_iszero(fx)) /* exact hit -> leave while-loop */ 144 break; 145 146 if (fx > 0.) { /* try another x */ 147 xtmp = x - step; 148 xtmp = _unur_max( xtmp, DISTR.domain[0] ); 149 } 150 else { 151 xtmp = x + step; 152 xtmp = _unur_min( xtmp, DISTR.domain[1] ); 153 } 154 155 fxtmp = CDF(xtmp) - U; 156 fxtmpabs = fabs(fxtmp); 157 158 if ( fxtmpabs < fxabs ) { /* improvement, update x */ 159 /* printf("fxabs: %g tmpabs: %g\n", fxabs, fxtmpabs); */ 160 step = 1.; /* set back stepsize */ 161 x = xtmp; 162 fx = fxtmp; 163 } 164 else if ( fxtmp*fx < 0. ) { /* step was too large, don't update x */ 165 step /= 2.; 166 } 167 else { /* step was too short, update x */ 168 step *= 2.; 169 x = xtmp; 170 fx = fxtmp; 171 } 172 173 dfx = PDF(x); 174 fxabs = fabs(fx); 175 176 if (flat_count < MAX_FLAT_COUNT) 177 flat_count++; 178 else { 179 _unur_error(gen->genid,UNUR_ERR_GEN_SAMPLING, 180 "Newton's method cannot leave flat region"); 181 x = _unur_max( x, DISTR.trunc[0]); 182 x = _unur_min( x, DISTR.trunc[1]); 183 return x; 184 } 185 } /* end of while-loop, (leaving flat region) */ 186 187 step = 1.; /* set back stepsize */ 188 189 if (_unur_iszero(fx)) /* exact hit -> finished */ 190 break; 191 192 193 if (_unur_isfinite(dfx)) { 194 do { /* newton-step (damped if nececcary) */ 195 damp /= 2.; 196 xtmp = x - damp * fx/dfx; 197 /* make sure that new point is inside (truncated) domain */ 198 xtmp = _unur_min( xtmp, DISTR.trunc[1] ); 199 xtmp = _unur_max( xtmp, DISTR.trunc[0] ); 200 fxtmp = CDF(xtmp) - U; 201 } while (fabs(fxtmp) > fxabs * (1.+UNUR_SQRT_DBL_EPSILON)); /* no improvement */ 202 } 203 else { 204 /* we cannot use Newton's rule if the derivative is not finite. */ 205 /* this happens when we hit a pole of the PDF. */ 206 /* use a bisection step instead. */ 207 xtmp = 0.5*(x + xold); 208 fxtmp = CDF(xtmp) - U; 209 } 210 211 /* updation variables according to newton-step */ 212 damp = 2.; /* set back factor for damping */ 213 xold = x; /* remember last value of x */ 214 x = xtmp; /* update x */ 215 fx = fxtmp; /* update function value at x */ 216 dfx = PDF(x); /* update derivative sof fx at x */ 217 fxabs = fabs(fx); /* update absolute value of fx */ 218 219 /* -- 2z. check stopping criterions -- */ 220 221 if ( GEN->x_resolution > 0. ) { 222 /* check x-error */ 223 /* we use a combination of absolute and relative x-error: */ 224 /* x-error < x-resolution * fabs(x) + x-resolution^2 */ 225 if ( _unur_iszero(fx) || /* exact hit */ 226 fabs(x-xold) < GEN->x_resolution * (fabs(x) + GEN->x_resolution) ) { 227 x_goal = TRUE; 228 } 229 else 230 x_goal = FALSE; 231 } 232 else { 233 /* no check */ 234 x_goal = TRUE; 235 } 236 237 if ( GEN->u_resolution > 0. ) { 238 /* check u-error */ 239 /* (we use a slightly smaller maximal tolerated error than given by user) */ 240 if ( fabs(fx) < 0.9 * rel_u_resolution ) { /* relative u resolution */ 241 u_goal = TRUE; 242 } 243 else if ( _unur_FP_same(xold, x) ) { 244 /* sharp peak or pole */ 245 _unur_warning(gen->genid,UNUR_ERR_GEN_SAMPLING, 246 "sharp peak or pole: accuracy goal in u cannot be reached"); 247 u_goal = TRUE; 248 } 249 else 250 u_goal = FALSE; 251 252 } 253 else { 254 u_goal = TRUE; 255 } 256 257 /* goal reached ? */ 258 if (x_goal && u_goal) 259 /*finished*/ 260 break; 261 } /* end of for-loop (MAXITER reached -> finished) */ 262 263 264 if (i >= GEN->max_iter) 265 _unur_warning(gen->genid,UNUR_ERR_GEN_SAMPLING, 266 "max number of iterations exceeded: accuracy goal might not be reached"); 267 268 /* make sure that result is within boundaries of (truncated) domain */ 269 x = _unur_max( x, DISTR.trunc[0]); 270 x = _unur_min( x, DISTR.trunc[1]); 271 272#ifdef UNUR_ENABLE_LOGGING 273 /* write info into LOG file (in case error) */ 274 if (gen->debug & NINV_DEBUG_SAMPLE) 275 _unur_ninv_debug_sample(gen, U, x, fx, i); 276#endif 277 278 return x; 279 280} /* end of _unur_ninv_sample_newton() */ 281 282/*---------------------------------------------------------------------------*/ 283 284#undef MAX_FLAT_COUNT 285 286/*---------------------------------------------------------------------------*/ 287