/* ------------------------------------------------------------------------- * * * Implementation of Algorithms HITRO-box and HITRO-plate * * Roman Karawatzki, Josef Leydold, and Klaus Poetzelberger (2006). * Automatic Markov Chain Monte Carlo Procedures for Sampling from * Multivariate Distributions. * * (c) 2006, Josef Leydold * * ------------------------------------------------------------------------- * * * ChangeLog: * * - Version 1.0 [Mon May 22 11:20:31] Josef Leydold: * First release * * ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------- * * * WARNING: These algorithms are Markov chain sampling algorithm, i.e. * they produce *correlated* random points. * * ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------- * * * DESCRIPTION: * * The HITRO (HIT-and-run + Ratio-Of-uniforms) Algorithm [ 'Hitro' is * the Slovenian word for fast ] combines the multivariate * Ratio-of-Uniforms method with the Hit-and-Run sampler. * * The Ratio-of-Uniforms transformation maps the region below the PDF, * i.e. * * G(PDF) = {(x,y): 0 < y < PDF(x)} * * into the region * * A(PDF) = {(u,v): 0 < v < (PDF(u/v^r + c))^(1/(r*n+1)) * * by means of the transformation * * R^n x (0,oo) -> R^n x (0,oo): * (u,v) -> (x,y) = (u/v^r + c, vr*n+1) * * where x = (x_1,...,x_n) and u = (u_1,...,u_n). * The vector 'c' should can be any point but should be chosen as a * "typical" point of the distribution (called 'center' below), * e.g. its mode. * The parameter 'r' controls the shape of the region A(PDF). * * For the algorithm it is important that the region A(PDF) is * bounded. The minimal bounding rectangle R(PDF) is given by * * v_min = 0 * v_max = sup PDF(x)^(1/(r*n+1)) * u_i,min = inf_{x_i} (x_i-c_i)*(PDF(x))^(r/(r*n+1)) * u_i,max = inf_{x_i} (x_i-c_i)*(PDF(x))^(r/(r*n+1)) * * Notice that A(PDF)is bounded if and only if * PDF(x) = O(||x||^(-n-1/r)). * * The value 'r'=1 is if special interest since then A(PDF) is convex * if and only if -1/PDF(x)^(-1/(n+1)) is concave. This is in * particular the case for log-concave densities. * * ------------------------------------------------------------------------- * * * APPLICATION PROGRAMMING INTERFACE: * * This is an implementation of the two algorithms HITRO-box and * HITRO-plate described in the above ACM TOMS paper but also adds * some of the variants. * * Arguments: * - number of variates (its "dimension") [ required ] * - probability density function (PDF) [ required ] * - center of the distribution [ required ] * ("typical" point of the distribution) * - parameter 'r' for Ratio-of-uniforms transformation [ required ] * (use '1.' if possible. for PDFs with higher tails * larger values are necessary) * - starting point for the chain [ optional ] * - uniform random number generator [ required ] * (the interface to this generator can be set by C macros) * * HITRO-plate: * - upper bound for the PDF [ optional ] * * HITRO-box: * - bounding rectangle for Ratio-of-Uniforms domain of [ optional ] * acceptance * - switch for 'adaptive uniform sampling' [ required ] * vs. 'simple rejection' * - switch for 'coordinate direction ("Gibbs") sampling' [ required ] * vs. 'random direction sampling' * * ------------------------------------------------------------------------- * * * HOW TO USE: * * (0) Use for favorite uniform random number generator. * Change type of and call to your uniform RNG in file * 'hitro.h' (if necessary). * * (1) Prepare a function of type 'double pdf(int dim, double* X)' * that evaluates the PDF at point 'X', where 'X' is a * double array of size 'dim'. The given function can return any * multiple of the PDF of the distribution. It is recommended, * that the PDF(mode) is about 1. * * (2) Decide whether to use HITRO-box or HITRO-plate. * HITRO-plate is simpler and faster when the bounding box is * not known. * * (3) Select 'r'. Set 'r'=1 if possible. * It is important that A(PDF) is bounded (the above infima * and supprema must be finite!) * If A(PDF) is not bounded use larger values for 'r'. * (But do not set 'r' much larger than necessary.) * * (4) Choose adaptive line sampling [Yes/No]. * Adaptive line sampling is much faster than simple rejection * samling. However, our proofs for convergence require 'r'=1 * and that the PDF is log-concave (see above). * There is strong experimental evidence that it also works for * non-logconcave (e.g. bimodal) densities as well (but it would * be better to run statistical tests on the output). * * Notice, that adaptive sampling is obligatory for HITRO-plate * and cannot be set by its interface. * * (5) Provide a "typical" point of the distribution, called 'center' * in the following, e.g. its mode. * Use NULL for the origin (0,...,0). * * (6) Provide an upper bound for the PDF (HITRO-plate) or the * left lower and right upper vertex of the bounding rectangle * (HITRO-box). Notice that the first coordinate is the * 'v'-coordinate, followed by the 'u' coordinates. * Furthermore, the bounding rectangle must be given for * PDF(X-'center'). * * These data are optional and can be omitted. Then an adaptive * approach is used to compute an upper bound and a bounding * rectangle ("slice sampling"). * Notice, however, that then two additional evaluations of the * PDF are necessary for each step of the chain. * * (7) Choose coordinate direction sampling ("Gibbs") or random * direction sampling. * We made the experience that coordinate direction sampling * performs better in Monte Carlo integration than random * direction sampling. * * Notice, that coordinate direction sampling is not available * for HITRO-plate. * * ------------------------------------------------------------------------- */ /* ------------------------------------------------------------------------- */ /* interface for uniform random number generator */ /* ------------------------------------------------------------------------- */ /* type of uniform random number generator */ typedef double (UNIFORM)(void); /* how to call uniform random number generator */ #define uniform(urng) urng(); /* ------------------------------------------------------------------------- */ /* API to HITRO Algorithms */ /* ------------------------------------------------------------------------- */ typedef struct hitro_chain *HITRO; /* HITRO chain */ /* Create chain for Algorithm HITRO-plate */ HITRO hitro_plate_create ( int dim, /* dimension of distribution. [ must be >= 1 ] */ double* center, /* center of distribution, typically the mode. [ array of size 'dim' ] [ use NULL for (0,...,0) ] */ double* start, /* starting point of chain. [ array of size 'dim' ] [ PDF(start) must not be > 0 ] [ use NULL if 'center'is used ] */ double pdfmax, /* upper bound for PDF. [ typically the PDF at the mode ] [ use 0. or negative value if not known ] */ double pdf(int dim, /* pointer to probability density function (PDF) */ double* X), /* [ must not return negative numbers ] [ 'X' is array of size 'dim' ] */ UNIFORM urng, /* pointer to uniform RNG */ double r /* parameter for Ratio-of-uniforms transformation */ ); /* Create chain for Algorithm HITRO-box */ HITRO hitro_box_create ( int dim, /* dimension of distribution. [ must be >= 1 ] */ double* center, /* center of distribution, typically the mode. [ array of size 'dim' ] [ use NULL for (0,...,0) ] */ double* start, /* starting point of chain. [ array of size 'dim' ] [ PDF(start) must not be > 0 ] [ use NULL if 'center'is used ] */ double* llbr, /* left lower point of bounding rectangle [ array of size 'dim'+1 ] [ first entry is 'v'-coordinate ] [ bounding rectangle for PDF(X-'center') ] [ use NULL if not known ] */ double* rubr, /* right upper point of bounding rectangle [ array of size 'dim'+1 ] [ first entry is 'v'-coordinate ] [ bounding rectangle for PDF(X-'center') ] [ use NULL if not known ] */ double pdf(int dim, /* pointer to probability density function (PDF) */ double* X), /* [ must not return negative numbers ] [ 'X' is array of size 'dim' ] */ UNIFORM urng, /* pointer to uniform RNG */ double r, /* parameter for Ratio-of-uniforms transformation */ int adaptive, /* 1 = use adaptive sampling from line segement 0 = use simple sampling */ int use_coord /* 1 = use coordinate direction ("Gibbs") sampling 0 = use random direction sampling */ ); /* Sample next point from HITRO chain */ int hitro_sample( HITRO chain, double *x ); /* Destroy HITRO chain */ int hitro_delete ( HITRO chain ); /* ------------------------------------------------------------------------- */