1 /***************************************************************************** 2 * * 3 * UNU.RAN -- Universal Non-Uniform Random number generator * 4 * * 5 ***************************************************************************** 6 * * 7 * FILE: pinv.c * 8 * * 9 * TYPE: continuous univariate random variate * 10 * METHOD: Polynomial interpolation based INVersion of CDF * 11 * * 12 * DESCRIPTION: * 13 * Compute values of CDF incrementally and interpolate resulting points * 14 * by polynomials. * 15 * * 16 * Integration: adaptive Gauss-Lobatto integration with 5 points * 17 * Interpolation: Newton recursion for interpolating polynomial * 18 * * 19 * REQUIRED: * 20 * pointer to the PDF, center of distribution * 21 * * 22 * OPTIONAL: * 23 * pointer to CDF, pointer to derivative of PDF * 24 * * 25 ***************************************************************************** 26 * * 27 * Copyright (c) 2008-2010 Wolfgang Hoermann and Josef Leydold * 28 * Department of Statistics and Mathematics, WU Wien, Austria * 29 * * 30 * This program is free software; you can redistribute it and/or modify * 31 * it under the terms of the GNU General Public License as published by * 32 * the Free Software Foundation; either version 2 of the License, or * 33 * (at your option) any later version. * 34 * * 35 * This program is distributed in the hope that it will be useful, * 36 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 37 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 38 * GNU General Public License for more details. * 39 * * 40 * You should have received a copy of the GNU General Public License * 41 * along with this program; if not, write to the * 42 * Free Software Foundation, Inc., * 43 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA * 44 * * 45 ***************************************************************************** 46 * * 47 * REFERENCES: * 48 * * 49 * [1] Derflinger, Gerhard, H{\"o}rmann, Wolfgang, and Leydold, Josef: * 50 * Random Variate Generation by Numerical Inversion when only the * 51 * Density Is Known. ACM Trans. Model. Comput. Simul. (2010), * 52 * to appear. * 53 * See also Research Report Series of the Department of * 54 * Statistics and Mathematics at WU Vienna, No. 90, June 2009. * 55 * http://epub.wu.ac.at/, oai:epub.wu-wien.ac.at:epub-wu-01_f41 * 56 * * 57 * [2] H{\"o}rmann, Wolfgang and Leydold, Josef: * 58 * Continuous Random Variate Generation by Fast Numerical Inversion, * 59 * ACM Trans. Model. Comput. Simul. 13(4), pp. 347--362, 2003. * 60 * * 61 ***************************************************************************** 62 * * 63 * Method PINV combines numerical integration with interpolation of the * 64 * inverse CDF. * 65 * * 66 * 1. Preprocessing: * 67 * * 68 * 1a. Estimate computationally relevant domain (support) of PDF * 69 * (finite interval where PDF is above some threshold value). * 70 * _unur_pinv_relevant_support() * 71 * _unur_pinv_searchborder() * 72 * * 73 * 1b. Compute area below PDF over relevant domain approximately. * 74 * _unur_pinv_approx_pdfarea() * 75 * * 76 * 1c. Compute computational domain where inverse CDF is approximated * 77 * (interval where we safely can compute coefficients of * 78 * interpolating polynomial). * 79 * _unur_pinv_computational_domain() * 80 * _unur_pinv_cut() * 81 * _unur_pinv_cut_bisect() * 82 * _unur_pinv_cut_CDF() * 83 * * 84 * 1d. Compute area below PDF over relevant domain with requested * 85 * accuracy and store subinterval boundaries and corresponding * 86 * from adaptive integration. * 87 * _unur_pinv_pdfarea() * 88 * * 89 * 2. Interpolation: * 90 * * 91 * 2a. Compute coefficients for interpolating polynomial for * 92 * fixed (sub-) interval. * 93 * _unur_pinv_newton_create() * 94 * _unur_pinv_chebyshev_points() * 95 * _unur_pinv_newton_cpoints() * 96 * * 97 * 2b. Evaluate interpolating polynomial. * 98 * _unur_pinv_newton_eval() * 99 * * 100 * 2c. Estimate approximation error for given interpolating polynomial. * 101 * _unur_pinv_newton_maxerror() * 102 * _unur_pinv_newton_testpoints() * 103 * _unur_pinv_cubic_hermite_monotone() * 104 * * 105 * * 106 * Currently the following methods are implemented: * 107 * * 108 * Quadrature (Integration): * 109 * Adaptive Gauss-Lobatto integration with 5 points. * 110 * see utils/lobatto.c * 111 * * 112 * Interpolation: * 113 * Newton recursion for coefficients of polynomial * 114 * ("Newton interpolation"). * 115 * * 116 * Hermite interpolation [2] (as the limit case with double and * 117 * tripple nodes). * 118 * * 119 *****************************************************************************/ 120 121 /*---------------------------------------------------------------------------*/ 122 123 /* enable additional code for development */ 124 /* #define PINV_DEVEL */ 125 126 /*---------------------------------------------------------------------------*/ 127 128 #include <unur_source.h> 129 #include <distr/distr.h> 130 #include <distr/distr_source.h> 131 #include <distr/cont.h> 132 #include <urng/urng.h> 133 #include <tests/unuran_tests.h> 134 #include <utils/lobatto_source.h> 135 #include "unur_methods_source.h" 136 #include "x_gen_source.h" 137 #include "pinv.h" 138 #include "pinv_struct.h" 139 140 141 /*---------------------------------------------------------------------------*/ 142 /* Constants */ 143 144 /* -- Global parameters */ 145 #define MAX_ORDER (17) 146 /* Maximum order of Newton interpolation polynomial */ 147 148 #define PINV_UERROR_CORRECTION (0.9) 149 /* PINV tries to create an approximation of the inverse CDF where the */ 150 /* U-error is bounded by the given u-resolution. However, the error caused */ 151 /* by cutting off the tails introduces some additional errors which must be */ 152 /* corrected. Thus the upper bound for the pure approximation error of the */ 153 /* interpolation is set to PINV_UERROR_CORRECTION * u-resolution. */ 154 155 #define PINV_DEFAULT_MAX_IVS (10000) 156 /* Default for maximum number of subintervals */ 157 158 /* -- Gauss-Lobatto integration */ 159 160 #define PINV_MAX_LOBATTO_IVS (20001) 161 /* Maximum number of subintervals for adaptive Gauss-Lobatto integration. */ 162 /* We keep this number fixed (independent of the maximum number of */ 163 /* subintervals for Newton interpolation), since the number of these */ 164 /* subintervals does neither depend on the order of the interpolating */ 165 /* polynomial nor on the requested u-resolution. */ 166 /* */ 167 /* Remark: This parameter MUST NOT be smaller than 2 !! */ 168 169 /* -- 1. Preprocessing */ 170 171 #define PINV_PDFLLIM (1.e-13) 172 /* Threshold value used for finding the boundary of the computational */ 173 /* domain. When starting the search at some point x0 then the search stops */ 174 /* when a point x is found where PDF(x) approx. PDF(x0) * PINV_PDFLLIM. */ 175 176 #define PINV_UERROR_AREA_APPROX (1.e-5) 177 /* Tolerated relative area when computing the area below the PDF */ 178 /* approximately in Step 1b. */ 179 180 #define PINV_TAILCUTOFF_FACTOR (0.05) 181 #define PINV_TAILCUTOFF_MAX (1.e-10) 182 /* For unbounded domains the tails has to be cut off. We use the given */ 183 /* u-resolution for finding the cut points. (The probability for each of the */ 184 /* chopped regions should be less than */ 185 /* HINV_TAILCUTOFF_FACTOR * u-resolution.) */ 186 /* However, the tail probabilities should not be greater than some threshold */ 187 /* value, given by PINV_TAILCUTOFF_MAX which reflects the precision of the */ 188 /* used stream of uniform pseudo-random numbers (typically about 2^32). */ 189 /* However, for computational reasons we use a value that is at least twice */ 190 /* the machine epsilon for the right hand boundary. */ 191 192 /* -- 2. Newton interpolation */ 193 194 #define PINV_UTOL_CORRECTION (0.05) 195 /* We use a smaller tolerance when computing the Gauss-Lobatto integral for */ 196 /* the PDF between construction points of the Newton polynomial. */ 197 198 #define PINV_MAX_ITER_IVS (10 * GEN->max_ivs) 199 /* maximum number of iterations for computing intervals for Newtwon */ 200 /* interpolation polynomials. Obviously it should be larger than the maximum */ 201 /* number of intervals (or the latter can be reduced). */ 202 203 #define PINV_GUIDE_FACTOR (1) 204 /* relative size of guide table for finding the subinterval corresponding */ 205 /* to the given U-value. */ 206 207 208 /*---------------------------------------------------------------------------*/ 209 /* Variants */ 210 211 #define PINV_VARIANT_PDF 0x0010u /* use PDF and Lobatto integration 212 [ if not present, use CDF ] */ 213 #define PINV_VARIANT_UPOINTS 0x0040u /* use Chebyshev points in u scale */ 214 #define PINV_VARIANT_KEEPCDF 0x0080u /* keep table for integration */ 215 216 /*---------------------------------------------------------------------------*/ 217 /* Debugging flags */ 218 /* bit 01 ... pameters and structure of generator (do not use here) */ 219 /* bits 02-12 ... setup */ 220 /* bits 13-24 ... adaptive steps */ 221 /* bits 25-32 ... trace sampling */ 222 223 #define PINV_DEBUG_REINIT 0x00000002u /* print parameters after reinit */ 224 #define PINV_DEBUG_TABLE 0x00000010u /* print table */ 225 #define PINV_DEBUG_SEARCHBD 0x00010000u /* trace search boundary */ 226 #define PINV_DEBUG_ITABLE 0x00020000u /* print table of integral values */ 227 228 /*---------------------------------------------------------------------------*/ 229 /* Flags for logging set calls */ 230 231 #define PINV_SET_ORDER 0x0001u /* order of polynomial */ 232 #define PINV_SET_ORDER_COR 0x1000u /* ... corrected */ 233 #define PINV_SET_SMOOTH 0x0002u /* smoothness of interpolant */ 234 #define PINV_SET_SMOOTH_COR 0x2000u /* ... corrected */ 235 #define PINV_SET_U_RESOLUTION 0x0004u /* maximal error in u */ 236 #define PINV_SET_UPOINTS 0x0008u /* use Chebyshev points in u scale */ 237 #define PINV_SET_BOUNDARY 0x0010u /* boundary of computational region */ 238 #define PINV_SET_SEARCHBOUNDARY 0x0020u /* search for boundary */ 239 #define PINV_SET_VARIANT 0x0040u /* variant of algorithm */ 240 #define PINV_SET_MAX_IVS 0x0080u /* maximum number of subintervals */ 241 #define PINV_SET_KEEPCDF 0x0100u /* keep table for integration */ 242 243 /*---------------------------------------------------------------------------*/ 244 245 #define GENTYPE "PINV" /* type of generator */ 246 247 /*---------------------------------------------------------------------------*/ 248 249 /*........................*/ 250 /* file: pinv_newset.ch */ 251 /*........................*/ 252 253 /* See pinv.h for 'new', 'set', and 'get' calls. */ 254 255 256 /*......................*/ 257 /* file: pinv_init.ch */ 258 /*......................*/ 259 260 static struct unur_gen *_unur_pinv_init (struct unur_par *par); 261 /*---------------------------------------------------------------------------*/ 262 /* Initialize new generator. */ 263 /*---------------------------------------------------------------------------*/ 264 265 /* static int _unur_pinv_reinit (struct unur_gen *gen); */ 266 /*---------------------------------------------------------------------------*/ 267 /* Reinitialize generator. */ 268 /*---------------------------------------------------------------------------*/ 269 270 static struct unur_gen *_unur_pinv_create (struct unur_par *par); 271 /*---------------------------------------------------------------------------*/ 272 /* create new (almost empty) generator object. */ 273 /*---------------------------------------------------------------------------*/ 274 275 static int _unur_pinv_check_par (struct unur_gen *gen); 276 /*---------------------------------------------------------------------------*/ 277 /* Check parameters of given distribution and method */ 278 /*---------------------------------------------------------------------------*/ 279 280 static struct unur_gen *_unur_pinv_clone (const struct unur_gen *gen); 281 /*---------------------------------------------------------------------------*/ 282 /* copy (clone) generator object. */ 283 /*---------------------------------------------------------------------------*/ 284 285 static void _unur_pinv_free (struct unur_gen *gen); 286 /*---------------------------------------------------------------------------*/ 287 /* destroy generator object. */ 288 /*---------------------------------------------------------------------------*/ 289 290 static int _unur_pinv_make_guide_table (struct unur_gen *gen); 291 /*---------------------------------------------------------------------------*/ 292 /* make a guide table for indexed search. */ 293 /*---------------------------------------------------------------------------*/ 294 295 static double _unur_pinv_eval_PDF (double x, struct unur_gen *gen); 296 /*---------------------------------------------------------------------------*/ 297 /* call to PDF. */ 298 /*---------------------------------------------------------------------------*/ 299 300 301 /*........................*/ 302 /* file: pinv_sample.ch */ 303 /*........................*/ 304 305 static double _unur_pinv_sample (struct unur_gen *gen); 306 /*---------------------------------------------------------------------------*/ 307 /* sample from generator */ 308 /*---------------------------------------------------------------------------*/ 309 310 static double _unur_pinv_eval_approxinvcdf (const struct unur_gen *gen, double u); 311 /*---------------------------------------------------------------------------*/ 312 /* evaluate interpolation of inverse CDF at u. */ 313 /*---------------------------------------------------------------------------*/ 314 315 /* declared in pinv.h: */ 316 /* double unur_pinv_eval_approxinvcdf (const struct unur_gen *gen, double u); */ 317 /* int unur_pinv_estimate_error (const UNUR_GEN *gen, int samplesize, double *max_error, double *MAE); */ 318 319 320 /*......................*/ 321 /* file: pinv_prep.ch */ 322 /*......................*/ 323 324 static int _unur_pinv_preprocessing (struct unur_gen *gen); 325 /*---------------------------------------------------------------------------*/ 326 /* 1. Find computational domain and compute PDF area. */ 327 /*---------------------------------------------------------------------------*/ 328 329 static int _unur_pinv_relevant_support (struct unur_gen *gen); 330 /*---------------------------------------------------------------------------*/ 331 /* 1a. Estimate computationally relevant domain (support) of PDF */ 332 /* (finite interval where PDF is above some threshold value). */ 333 /*---------------------------------------------------------------------------*/ 334 335 static double _unur_pinv_searchborder (struct unur_gen *gen, double x0, double bound, 336 double *dom, int *search); 337 /*---------------------------------------------------------------------------*/ 338 /* [1a.] find left or right hand border of relevant domain. */ 339 /*---------------------------------------------------------------------------*/ 340 341 static int _unur_pinv_approx_pdfarea (struct unur_gen *gen); 342 /*---------------------------------------------------------------------------*/ 343 /* 1b. Compute area below PDF over relevant domain approximately. */ 344 /*---------------------------------------------------------------------------*/ 345 346 static int _unur_pinv_pdfarea (struct unur_gen *gen); 347 /*---------------------------------------------------------------------------*/ 348 /* 1d. Compute area below PDF with requested accuracy and */ 349 /* store intermediate results from adaptive integration. */ 350 /*---------------------------------------------------------------------------*/ 351 352 static int _unur_pinv_computational_domain (struct unur_gen *gen); 353 /*---------------------------------------------------------------------------*/ 354 /* 1c. Compute computational domain where inverse CDF is approximated */ 355 /* (interval where we safely can compute coefficients of */ 356 /* interpolating polynomial). */ 357 /*---------------------------------------------------------------------------*/ 358 359 static double _unur_pinv_cut (struct unur_gen *gen, double w, double dw, double crit); 360 /*---------------------------------------------------------------------------*/ 361 /* [1c.] calculate cut-off points for computational domain of distribution. */ 362 /*---------------------------------------------------------------------------*/ 363 364 static double _unur_pinv_cut_bisect (struct unur_gen *gen, double x0, double x1); 365 /*---------------------------------------------------------------------------*/ 366 /* [1c.] calculate cut-off points as boudary of bounded support of PDF. */ 367 /*---------------------------------------------------------------------------*/ 368 369 static int _unur_pinv_computational_domain_CDF (struct unur_gen *gen); 370 /*---------------------------------------------------------------------------*/ 371 /* 1c. Compute computational domain where inverse CDF is approximated */ 372 /* (interval where we safely can compute coefficients of */ 373 /* interpolating polynomial). */ 374 /* Use CDF. */ 375 /*---------------------------------------------------------------------------*/ 376 377 static double _unur_pinv_cut_CDF( struct unur_gen *gen, double dom, double x0, double ul, double uu ); 378 /*---------------------------------------------------------------------------*/ 379 /* [1c.] calculate cut-off points for computational domain of distribution */ 380 /* using CDF. */ 381 /*---------------------------------------------------------------------------*/ 382 383 static double _unur_pinv_Udiff (struct unur_gen *gen, double x, double h, double *fx); 384 /*---------------------------------------------------------------------------*/ 385 /* compute difference CDF(x+h)-CDF(x) (approximately), where CDF is the */ 386 /* integral of the given (quasi-) density. */ 387 /*---------------------------------------------------------------------------*/ 388 389 390 /*........................*/ 391 /* file: pinv_newton.ch */ 392 /*........................*/ 393 394 static int _unur_pinv_create_table( struct unur_gen *gen ); 395 /*---------------------------------------------------------------------------*/ 396 /* create table for Newton interpolation */ 397 /*---------------------------------------------------------------------------*/ 398 399 static int _unur_pinv_chebyshev_points (double *pt, int order, int smooth); 400 /*---------------------------------------------------------------------------*/ 401 /* [2a.] Compute Chebyshev points. */ 402 /*---------------------------------------------------------------------------*/ 403 404 static int _unur_pinv_newton_cpoints (double *xval, int order, struct unur_pinv_interval *iv, 405 double h, double *chebyshev, int smooth, int use_upoints); 406 /*---------------------------------------------------------------------------*/ 407 /* [2a.] Compute points for construct interpolating polynomial. */ 408 /*---------------------------------------------------------------------------*/ 409 410 static int _unur_pinv_newton_create (struct unur_gen *gen, struct unur_pinv_interval *iv, 411 double *xval); 412 /*---------------------------------------------------------------------------*/ 413 /* 2a. Compute coefficients for Newton interpolation. */ 414 /*---------------------------------------------------------------------------*/ 415 416 static int _unur_pinv_linear_create (struct unur_gen *gen, struct unur_pinv_interval *iv, 417 double *xval); 418 /*---------------------------------------------------------------------------*/ 419 /* [2a.] Compute coefficients for linear interpolation. */ 420 /*---------------------------------------------------------------------------*/ 421 422 static double _unur_pinv_newton_eval (double q, double *ui, double *zi, int order); 423 /*---------------------------------------------------------------------------*/ 424 /* 2b. evaluate Newton polynomial. */ 425 /*---------------------------------------------------------------------------*/ 426 427 static double _unur_pinv_newton_maxerror (struct unur_gen *gen, struct unur_pinv_interval *iv, 428 double *xval, int use_linear); 429 /*---------------------------------------------------------------------------*/ 430 /* 2c. estimate maximal error of Newton interpolation in subinterval */ 431 /*---------------------------------------------------------------------------*/ 432 433 static int _unur_pinv_newton_testpoints (double *utest, double ui[], int order); 434 /*---------------------------------------------------------------------------*/ 435 /* [2c.] calculate the local maxima of the interpolation polynomial */ 436 /*---------------------------------------------------------------------------*/ 437 438 static int _unur_pinv_linear_testpoints (double *utest, double *ui, int order); 439 /*---------------------------------------------------------------------------*/ 440 /* [2c.] create table of test points for linear interpolation */ 441 /*---------------------------------------------------------------------------*/ 442 443 static int _unur_pinv_cubic_hermite_is_monotone(); 444 /*---------------------------------------------------------------------------*/ 445 /* [2c.] check monotonicity of cubic Hermite interpolation */ 446 /*---------------------------------------------------------------------------*/ 447 448 static int _unur_pinv_interval( struct unur_gen *gen, int i, double x, double cdfx ); 449 /*---------------------------------------------------------------------------*/ 450 /* make a new interval i with left boundary point x and CDF(x). */ 451 /*---------------------------------------------------------------------------*/ 452 453 static int _unur_pinv_lastinterval( struct unur_gen *gen ); 454 /*---------------------------------------------------------------------------*/ 455 /* update size of array and set all uninitialized values to 0. */ 456 /*---------------------------------------------------------------------------*/ 457 458 459 /*.......................*/ 460 /* file: pinv_debug.ch */ 461 /*.......................*/ 462 463 #ifdef UNUR_ENABLE_LOGGING 464 /*---------------------------------------------------------------------------*/ 465 /* the following functions print debugging information on output stream, */ 466 /* i.e., into the LOG file if not specified otherwise. */ 467 /*---------------------------------------------------------------------------*/ 468 469 static void _unur_pinv_debug_init_start (const struct unur_gen *gen); 470 /*---------------------------------------------------------------------------*/ 471 /* print before setup starts. */ 472 /*---------------------------------------------------------------------------*/ 473 474 static void _unur_pinv_debug_init (const struct unur_gen *gen, int ok); 475 /*---------------------------------------------------------------------------*/ 476 /* print after generator has been initialized has completed. */ 477 /*---------------------------------------------------------------------------*/ 478 479 static void _unur_pinv_debug_relevant_support (const struct unur_gen *gen); 480 /*---------------------------------------------------------------------------*/ 481 /* print relevant domain */ 482 /*---------------------------------------------------------------------------*/ 483 484 static void _unur_pinv_debug_pdfarea (const struct unur_gen *gen, int approx); 485 /*---------------------------------------------------------------------------*/ 486 /* print estimated area below PDF */ 487 /*---------------------------------------------------------------------------*/ 488 489 static void _unur_pinv_debug_computational_domain (const struct unur_gen *gen); 490 /*---------------------------------------------------------------------------*/ 491 /* print computational domain */ 492 /*---------------------------------------------------------------------------*/ 493 494 static void _unur_pinv_debug_intervals (const struct unur_gen *gen ); 495 /*---------------------------------------------------------------------------*/ 496 /* print starting points or table for algorithms into LOG file. */ 497 /*---------------------------------------------------------------------------*/ 498 499 static void _unur_pinv_debug_create_table (const struct unur_gen *gen, 500 int iter, int n_incr_h, int n_decr_h, 501 int n_use_linear); 502 /*---------------------------------------------------------------------------*/ 503 /* print data that have been collected while creating polynomials. */ 504 /*---------------------------------------------------------------------------*/ 505 #endif 506 507 508 /*......................*/ 509 /* file: pinv_info.ch */ 510 /*......................*/ 511 512 #ifdef UNUR_ENABLE_INFO 513 static void _unur_pinv_info( struct unur_gen *gen, int help ); 514 /*---------------------------------------------------------------------------*/ 515 /* create info string. */ 516 /*---------------------------------------------------------------------------*/ 517 #endif 518 519 /*---------------------------------------------------------------------------*/ 520 /* abbreviations */ 521 522 #define DISTR_IN distr->data.cont /* data for distribution object */ 523 524 #define PAR ((struct unur_pinv_par*)par->datap) /* data for parameter object */ 525 #define GEN ((struct unur_pinv_gen*)gen->datap) /* data for generator object */ 526 #define DISTR gen->distr->data.cont /* data for distribution in generator object */ 527 528 #define SAMPLE gen->sample.cont /* pointer to sampling routine */ 529 530 #define PDF(x) (_unur_pinv_eval_PDF((x),(gen))) /* call to PDF */ 531 #define dPDF(x) (_unur_cont_dPDF((x),(gen->distr))) /* call to derivate of PDF */ 532 #define CDF(x) (_unur_cont_CDF((x),(gen->distr))) /* call to CDF */ 533 534 /*---------------------------------------------------------------------------*/ 535 536 #define _unur_pinv_getSAMPLE(gen) (_unur_pinv_sample) 537 538 /*---------------------------------------------------------------------------*/ 539 540 541 542 /*---------------------------------------------------------------------------*/ 543 /* since there is only file scope or program code, we abuse the */ 544 /* #include directive. */ 545 546 /** Public: User Interface (API) **/ 547 #include "pinv_newset.ch" 548 549 /** Private **/ 550 #include "pinv_init.ch" 551 #include "pinv_sample.ch" 552 #include "pinv_prep.ch" 553 #include "pinv_newton.ch" 554 #include "pinv_debug.ch" 555 #include "pinv_info.ch" 556 557 /*---------------------------------------------------------------------------*/ 558