1This is ranlip.info, produced by makeinfo version 4.5 from 2ranlip.texinfo. 3 4 =insert GNU GENERAL PUBLIC LICENSE Version 2, 5June 1991 6 7 Copyright (C) 1989, 1991 Free Software Foundation, Inc. 59 Temple 8Place, Suite 330, Boston, MA 02111-1307 USA. Everyone is permitted to 9copy and distribute verbatim copies of this license document, but 10changing it is not allowed. 11 12 13 14File: ranlip.info, Node: Top, Next: Introduction, Up: (dir) 15 16Class Library ranlip for Multivariate Non-uniform Random Variate Generation 17*************************************************************************** 18 19 This manual describes generation of non-uniform random variates from 20Lipschitz-continuous densities using acceptance/ rejection, and the 21class library *ranlip* which implements this method. It is assumed that 22the required distribution has Lipschitz-continuous density, which is 23either given analytically or as a black box. The algorithm builds a 24piecewise constant upper approximation to the density (the hat 25function), using a large number of its values and subdivision of the 26domain into hyper rectangles. 27 28 The class library *ranlip* provides very competitive preprocessing 29and generation times, and yields small rejection constant, which is a 30measure of efficiency of the generation step. It exhibits good 31performance for up to five variables, and provides the user with a 32black box non-uniform random variate generator for a large class of 33distributions, in particular, multimodal distributions. 34 35 GNU GENERAL PUBLIC LICENSE Version 2, 36June 1991 37 38 Copyright (C) 1989, 1991 Free Software Foundation, Inc. 59 Temple 39Place, Suite 330, Boston, MA 02111-1307 USA. Everyone is permitted to 40copy and distribute verbatim copies of this license document, but 41changing it is not allowed. 42 The menu below lists the major sections which give a brief 43background, overview of the library and illustrative examples on how to 44use it. 45 46* Menu: 47 48* Introduction:: Overview of software. 49* Description of Library:: Library interface methods. 50* Examples:: Library in action. 51* Performance:: Computational complexity and performance of the algorithms. 52* Index:: Complete index. 53 54 55File: ranlip.info, Node: Introduction, Next: Description of Library, Prev: Top, Up: Top 56 57Introduction 58************ 59 60 This manual describes the programming library *ranlip*, which 61implements the method of acceptance/ rejection in the multivariate 62case, for Lipschitz continuous densities. It assumes that the 63Lipschitz constant of the density rho is known, or can be approximated, 64and that computation of the values of rho at distinct points is not 65expensive. The method builds a piecewise constant hat function, by 66subdividing the domain into hyper rectangles, and by using a large 67number of values of rho. Lipschitz properties of rhoallow one to 68overestimate rho at all other points, and thus to overestimate the 69absolute maxima of rho on the elements of the partition. 70 71 The class library *ranlip* implements computation of the hat 72function and generation of random variates, and makes this process 73transparent to the user. The user needs to provide a method of 74evaluation of rho at a given point, and the number of elements in the 75subdivision of the domain, which is the parameter characterizing the 76quality of the hat function and the number of computations at the 77preprocessing step. 78 79 The class of Lipschitz-continuous densities is very broad, and 80includes many multimodal densities, which are hard to deal with. No 81other properties beyond Lipschitz continuity are required, and the 82Lipschitz constant, if not provided, can be estimated automatically. 83The algorithm does not require rho to be given analytically, to be 84differentiable, or to be normalized. 85 86 87File: ranlip.info, Node: Description of Library, Next: Examples, Prev: Introduction, Up: Top 88 89Description of Library 90********************** 91 92 The main class which provides the interface to the preprocessing and 93computation is called `CRanLip'. It is illustrated in the following 94sections, together with amore extensive description of the interface. 95Can be viewed via the following menu: 96 97* Menu: 98 99* Class CRanLip:: interface of class CRanLip, inheriting from CRanLip 100* Closer Look at Interface:: more detailed view of the interface. 101 102 103File: ranlip.info, Node: Class CRanLip, Next: Closer Look at Interface, Up: Description of Library 104 105Class CRanLip 106============= 107 108 The main class which provides the interface to the preprocessing and 109random variate generation is called CRanLip. This is an abstract class, 110from which the user must derive his own class which overrides the 111method for rho(x), and declare an instance of that class. The interface 112of Class CRandLip is illustrated bellow followed by an example of how 113to inherit from CRanLip. 114 115 class CRanLip { 116 public: 117 // initialises the tables and arrays for the generator 118 // should be the first method to call 119 void Init(int dim, double* left, double* right); 120 121 // sets the pointer to the uniform random number generator. 122 // The default is ranlux generator 123 void SetUniformGenerator(UFunction gen); 124 125 // sets the seed for the uniform random number generator 126 void Seed(int seed); 127 128 // computes the hat function given the Lipschitz constant 129 // num and numfine determine the rough and fine partitions 130 void PrepareHatFunction(int num, int numfine, double Lip); 131 132 // computes the hat function, and automatically computes 133 // the Lipschitz constant 134 void PrepareHatFunctionAuto(int num, int numfine, double minLip); 135 136 // generates a random variate with the required density 137 void RandomVec(double* p); 138 139 // frees the memory occupied by the partition and various tables 140 void FreeMem(); 141 ... 142 } 143 144 The example bellow illustrates how the user needs to inherit from 145CRanLip. 146 147 // declare a derived class MyRnumGen, with one method to override 148 class MyRnumGen:public CRanLip { 149 public: virtual double Distribution(double* p) ; 150 }; 151 152 double MyRnumGen::Distribution(double* p) 153 { // example: multivriate normal distribution 154 double r; 155 for(int j=0;j<Dimension;j++) { 156 r+=p[j]*p[j]; 157 } 158 return exp(-r);; 159 } 160 // declare an instance of this class 161 MyRnumGen MyGen; 162 163 164File: ranlip.info, Node: Closer Look at Interface, Prev: Class CRanLip, Up: Description of Library 165 166Closer Look at Interface 167======================== 168 169Init(int dim, double* left, double* right) 170------------------------------------------ 171 172 Initialises the internal variables of this class. dim is the 173dimension, left and right are arrays of size dim which determine the 174domain of rho: left_i <= x_i <= right_i. *Init must be called only once 175before any other method.* 176 177SetUniformGenerator(UFunction gen) 178---------------------------------- 179 180 Sets a pointer to the uniform random number generator on (0,1). The 181default is M. Luescher's ranlux generator, but the user can override it 182and use his own preferred generator. 183 184Seed(int seed) 185-------------- 186 187 Sets the seed of the default uniform random number generator ranlux. 188If the user has supplied his own generator in SetUniformGenerator, that 189generator's seed() function should be called instead. 190 191PrepareHatFunction(int num, int numfine, double Lip) 192---------------------------------------------------- 193 194 Builds the hat function, using the Lipschitz constant supplied in 195Lip. Parameters num and numfine determine the rough and fine 196partitions. num is the number of subdivisions in each variable to 197partition the domain D into hyper rectangles D_k. On each D_k, the hat 198function will have a constant value h_k. numfine(>1) is the number of 199subdivisions in the finer partition in each variable. Each set D_k is 200subdivided into (numfine-1)^dim smaller hyper rectangles, in order to 201improve the quality of the overestimate h_k. There will be in total 202(num*numfine)^dim evaluations of rho (calls to Distribution()). numfine 203should be a power of 2 for numerical efficiency reasons (if not, it 204will be automatically changed to a power of 2 larger than the supplied 205value). numfine can be 2, in which case the fine partition is not used. 206 207PrepareHatFunctionAuto(int num, int numfine, double minLip=0) 208------------------------------------------------------------- 209 210 Builds the hat function and automatically computes an estimate to 211the Lipschitz constant. Parameters num and numfine determine the rough 212and fine partitions, and are described in PrepareHatFunction() method. 213minLip denotes the lower bound on the value of the computed Lipschitz 214constant, the default value is 0. 215 216RandomVec(double* p) 217-------------------- 218 219 Generates a random variate with the density rho. Should be called 220after PrepareHatFunction() or PrepareHatFunctionAuto(). The parameter p 221is an array of size dim, it will contain the components of the computed 222random vector. 223 224FreeMemory() 225------------ 226 227 Frees the memory occupied by the data structures, which can be very 228large. It destroys the hat function, and RandomVec() method cannot be 229called after FreeMemory(). Automatically called from the destructor. 230This method is useful to deallocate memory while the object CRanLip 231still exists. 232 233 234File: ranlip.info, Node: Examples, Next: Performance, Prev: Description of Library, Up: Top 235 236Exemples 237******** 238 239 There are several examples of the usage of *ranlip* provided in the 240distribution. There are three basic steps: to supply the required 241distribution, to build the hat function, and to generate random 242variates. As illustrated in following examples. 243 244* Menu: 245 246* Common Use Example:: Common usage. 247* Another Example:: Common usage. 248* Procedural Example:: Uses C syntax. 249 250 251File: ranlip.info, Node: Common Use Example, Next: Another Example, Up: Examples 252 253Common Use Example 254================== 255 256 #include "ranlip.h" 257 #define dim 4 // the dimension 258 // declare a derived class MyRnumGen, with one method to override 259 class MyRnumGen:public CRanLip { 260 public: virtual double Distribution(double* p) ; 261 }; 262 double MyRnumGen::Distribution(double* p) 263 { // example: multivriate normal distribution 264 double r=0.0; 265 for(int j=0;j<Dimension;j++) r+=p[j]*p[j]; 266 return exp(-r); 267 } 268 void main(int argc, char *argv[]){ 269 double LipConst = 4.0; 270 MyRnumGen MyGen; 271 double left[dim], right[dim], p[dim]; 272 int i; 273 // set the domain to be the unit hypercube 274 for(i=0;i<dim;i++) {left[i]=0; right[i]=1;} 275 276 MyGen.Init(dim,left,right); 277 MyGen.PrepareHatFunction(10,8,LipConst); 278 MyGen.Seed(10); 279 for(i=0;i<1000;i++) { 280 MyGen.RandomVec(p); 281 // do something with p 282 } 283 MyGen.FreeMemory(); 284 // now MyGen can be reused 285 MyGen.Init(dim,left,right); 286 MyGen.PrepareHatFunctionAuto(10,32); 287 for(i=0;i<1000;i++) { 288 MyGen.RandomVec(p); 289 } 290 } 291 292 293File: ranlip.info, Node: Another Example, Next: Procedural Example, Prev: Common Use Example, Up: Examples 294 295Another Example 296=============== 297 298 // Example 2 299 #include "ranlip.h" 300 #define dim 3 // the dimension 301 // declare a derived class MyRnumGen, with one method to override 302 class MyRnumGen:public CRanLip { 303 public: virtual double Distribution(double* p) ; 304 }; 305 double MyRnumGen::Distribution(double* p) 306 { // example: multivriate normal distribution 307 double r=0.0; 308 for(int j=0;j<Dimension;j++) r+=p[j]*p[j]; 309 return exp(-r); 310 } 311 void main(int argc, char *argv[]){ 312 double LipConst; 313 MyRnumGen MyGen; 314 double left[dim], right[dim], p[dim]; 315 int i; 316 // set the domain to be a hypercube 317 for(i=0;i<dim;i++) {left[i]=-2; right[i]=2;} 318 MyGen.Init(dim,left,right); 319 MyGen.PrepareHatFunctionAuto(10,32,0.01); 320 MyGen.Seed(10); 321 322 for(i=0;i<1000;i++) { 323 MyGen.RandomVec(p); 324 // do something with p 325 } 326 cout<<"acceptance ratio is "<<1000.0/ MyGen.count_total<<endl; 327 LipConst=MyGen.Lipschitz; // computed Lipschitz constant 328 if(MyGen.count_error>0) { // Lipschitz constant was too low 329 LipConst*=2; 330 MyGen.PrepareHatFunction(10,32,LipConst); 331 } 332 for(i=0;i<1000;i++) { 333 MyGen.RandomVec(p); 334 // do something with p 335 } 336 } 337 338 339File: ranlip.info, Node: Procedural Example, Prev: Another Example, Up: Examples 340 341Procedural Example 342================== 343 344 // Example 3 using procedural interface 345 #include "ranlipproc.h" 346 #define Dim 3 // the dimension 347 // implement calculation of the density in MyDist function 348 double MyDist(double* p, int dim) 349 { // example: multivariate normal distribution 350 double r=0.0; 351 for(int j=0;j<dim;j++) r+=p[j]*p[j]; 352 return exp(-r); 353 } 354 double MyRand() // use my own random number generator 355 { return (double)rand()/(RAND_MAX+0.001); } //random numbers in [0,1) 356 357 void main(int argc, char *argv[]){ 358 double LipConst; int i; 359 double left[Dim], right[Dim], p[Dim]; 360 // set the domain to be a hypercube 361 for(i=0;i<Dim;i++) {left[i]=-2; right[i]=2;} 362 363 InitRanLip(Dim,a,b); 364 // pass the address of the density function (required) 365 SetDistFunctionRanLip(&MyDist); 366 // pass the address of my generator (optional) 367 SetUniformGeneratorRanLip(&MyRand); 368 srand(10); // use srand, not SeedRanLip 369 370 PrepareHatFunctionAutoRanLip(10,8); 371 for(i=0;i<1000;i++) { 372 RandomVecRanLip(p); 373 // do something with p 374 } 375 cout<<"acceptance ratio is "<<1000.0/ Count_totalRanLip() <<endl; 376 LipConst=LipschitzRanLip(); // computed Lipschitz constant 377 if(Count_errorRanLip()>0) { // Lipschitz constant was too low 378 LipConst*=2; 379 PrepareHatFunctionRanLip(10,32,LipConst); 380 } 381 for(i=0;i<1000;i++) 382 RandomVecRanLip(p); 383 } 384 385 386File: ranlip.info, Node: Performance, Next: Index, Prev: Examples, Up: Top 387 388Performance 389*********** 390 391 It is important to estimate the computing time and memory 392requirements when using *ranlip*, especially for the case of several 393variables. The quality of the hat function directly depends on the 394number of values of rho(x) used for its computation. The higher this 395number is, the longer is preprocessing step (building the hat 396function), but the more efficient is the generation step (less rejected 397variates). A good estimate of the Lipschitz constant is also important, 398as it improves the quality of the hat function. 399 400 The method PrepareHatFunction(num, numfine, LipConst) uses the first 401two parameters to establish the rough partition of the domain, sets 402D_k, on which the hat function will have a constant value h_k, and the 403fine partition, used to compute this value. In total, (num*numfine)^dim 404evaluations of rho will be performed. The memory required by the 405algorihm is numfine^dim + 2*num^dim values of type double (8 bytes 406each). 407 408 For numerical efficiency, the second parameter numfine should be a 409power of 2. This facilitates computation of the neighbours in an 410n-dimensional mesh by using binary arithmetic. 411 412 413File: ranlip.info, Node: Index, Prev: Performance, Up: Top 414 415Index 416***** 417 418* Menu: 419 420* chapter, Computational complexity and performance: Performance. 421* chapter, exemples: Examples. 422* chapter, Introduction: Introduction. 423* chapter, ranlip Description: Description of Library. 424* Examples, Another Example: Another Example. 425* Examples, Common Use Example: Common Use Example. 426* Examples, Procedural Example: Procedural Example. 427* section, interface of class CRanLip: Class CRanLip. 428* section, RanLip library Interface: Closer Look at Interface. 429 430 431 432Tag Table: 433Node: Top388 434Node: Introduction2228 435Node: Description of Library3805 436Node: Class CRanLip4348 437Node: Closer Look at Interface6457 438Node: Examples9384 439Node: Common Use Example9871 440Node: Another Example11140 441Node: Procedural Example12618 442Node: Performance14270 443Node: Index15507 444 445End Tag Table 446