1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: gibbs.c *
8 * *
9 * TYPE: continuous multivariate random variate *
10 * METHOD: Gibbs sampler using full conditional distributions. *
11 * *
12 * DESCRIPTION: *
13 * Given PDF *
14 * Produce a value x consistent with its density *
15 * *
16 * REQUIRED: *
17 * pointer to the density function and its derivatives *
18 * *
19 * OPTIONAL: *
20 * mode of the density *
21 * *
22 *****************************************************************************
23 * *
24 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
25 * Department of Statistics and Mathematics, WU Wien, Austria *
26 * *
27 * This program is free software; you can redistribute it and/or modify *
28 * it under the terms of the GNU General Public License as published by *
29 * the Free Software Foundation; either version 2 of the License, or *
30 * (at your option) any later version. *
31 * *
32 * This program is distributed in the hope that it will be useful, *
33 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
34 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
35 * GNU General Public License for more details. *
36 * *
37 * You should have received a copy of the GNU General Public License *
38 * along with this program; if not, write to the *
39 * Free Software Foundation, Inc., *
40 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
41 * *
42 *****************************************************************************
43 * *
44 * REFERENCES: *
45 * *
46 *****************************************************************************/
47
48 /*---------------------------------------------------------------------------*/
49
50 #include <unur_source.h>
51 #include <distr/distr.h>
52 #include <distr/distr_source.h>
53 #include <distributions/unur_distributions.h>
54 #include <distr/condi.h>
55 #include <distr/cvec.h>
56 #include <urng/urng.h>
57 #include <utils/matrix_source.h>
58 #include "unur_methods_source.h"
59 #include "x_gen.h"
60 #include "x_gen_source.h"
61 #include "arou.h"
62 #include "ars.h"
63 #include "tdr.h"
64
65 #include "gibbs.h"
66 #include "gibbs_struct.h"
67
68 #ifdef UNUR_ENABLE_INFO
69 # include <tests/unuran_tests.h>
70 #endif
71
72 /*---------------------------------------------------------------------------*/
73 /* Variants */
74
75 #define GIBBS_VARMASK_VARIANT 0x000fu /* indicates variant */
76 #define GIBBS_VARIANT_COORD 0x0001u /* coordinate sampler */
77 #define GIBBS_VARIANT_RANDOMDIR 0x0002u /* random direction sampler */
78
79 #define GIBBS_VARMASK_T 0x00f0u /* indicates transformation */
80 #define GIBBS_VAR_T_SQRT 0x0010u /* T(x) = -1/sqrt(x) */
81 #define GIBBS_VAR_T_LOG 0x0020u /* T(x) = log(x) */
82 #define GIBBS_VAR_T_POW 0x0030u /* T(x) = -x^c */
83
84 /*---------------------------------------------------------------------------*/
85 /* Debugging flags */
86 /* bit 01 ... pameters and structure of generator (do not use here) */
87 /* bits 02-12 ... setup */
88 /* bits 13-24 ... adaptive steps */
89 /* bits 25-32 ... trace sampling */
90
91 #define GIBBS_DEBUG_CONDI 0x01000000u
92
93 /*---------------------------------------------------------------------------*/
94 /* Flags for logging set calls */
95
96 #define GIBBS_SET_C 0x001u /* set parameter for transformation T */
97 #define GIBBS_SET_X0 0x002u /* set starting point */
98 #define GIBBS_SET_THINNING 0x004u /* set thinning factor */
99 #define GIBBS_SET_BURNIN 0x008u /* set length of burn-in */
100
101 /*---------------------------------------------------------------------------*/
102
103 #define GENTYPE "GIBBS" /* type of generator */
104
105 /*---------------------------------------------------------------------------*/
106
107 static struct unur_gen *_unur_gibbs_init( struct unur_par *par );
108 static int _unur_gibbs_coord_init( struct unur_gen *gen );
109 static int _unur_gibbs_randomdir_init( struct unur_gen *gen );
110 /*---------------------------------------------------------------------------*/
111 /* Initialize new generator. */
112 /*---------------------------------------------------------------------------*/
113
114 static struct unur_gen *_unur_gibbs_create( struct unur_par *par );
115 /*---------------------------------------------------------------------------*/
116 /* create new (almost empty) generator object. */
117 /*---------------------------------------------------------------------------*/
118
119 static struct unur_gen *_unur_gibbs_clone( const struct unur_gen *gen );
120 /*---------------------------------------------------------------------------*/
121 /* copy (clone) generator object. */
122 /*---------------------------------------------------------------------------*/
123
124 static void _unur_gibbs_free( struct unur_gen *gen);
125 /*---------------------------------------------------------------------------*/
126 /* destroy generator object. */
127 /*---------------------------------------------------------------------------*/
128
129 static int _unur_gibbs_coord_sample_cvec( struct unur_gen *gen, double *vec );
130 static int _unur_gibbs_randomdir_sample_cvec( struct unur_gen *gen, double *vec );
131 /*---------------------------------------------------------------------------*/
132 /* sample from generator */
133 /*---------------------------------------------------------------------------*/
134
135 static struct unur_gen *_unur_gibbs_normalgen( struct unur_gen *gen );
136 /*---------------------------------------------------------------------------*/
137 /* create a normal random variate generator */
138 /*---------------------------------------------------------------------------*/
139
140 static void _unur_gibbs_random_unitvector( struct unur_gen *gen, double *direction );
141 /*---------------------------------------------------------------------------*/
142 /* generate a random direction vector */
143 /*---------------------------------------------------------------------------*/
144
145 #ifdef UNUR_ENABLE_LOGGING
146 /*---------------------------------------------------------------------------*/
147 /* the following functions print debugging information on output stream, */
148 /* i.e., into the LOG file if not specified otherwise. */
149 /*---------------------------------------------------------------------------*/
150
151 static void _unur_gibbs_debug_init_start( const struct unur_gen *gen );
152 /*---------------------------------------------------------------------------*/
153 /* print before init of generator starts. */
154 /*---------------------------------------------------------------------------*/
155
156 static void _unur_gibbs_debug_init_condi( const struct unur_gen *gen );
157 /*---------------------------------------------------------------------------*/
158 /* print list of conditional generators. */
159 /*---------------------------------------------------------------------------*/
160
161 static void _unur_gibbs_debug_burnin_failed( const struct unur_gen *gen );
162 /*---------------------------------------------------------------------------*/
163 /* print after burnin has failed. */
164 /*---------------------------------------------------------------------------*/
165
166 static void _unur_gibbs_debug_init_finished( const struct unur_gen *gen, int success );
167 /*---------------------------------------------------------------------------*/
168 /* print after generator has been initialized has completed. */
169 /*---------------------------------------------------------------------------*/
170 #endif
171
172 #ifdef UNUR_ENABLE_INFO
173 static void _unur_gibbs_info( struct unur_gen *gen, int help );
174 /*---------------------------------------------------------------------------*/
175 /* create info string. */
176 /*---------------------------------------------------------------------------*/
177 #endif
178
179 /*---------------------------------------------------------------------------*/
180 /* abbreviations */
181
182 #define DISTR_IN distr->data.cvec /* data for distribution object */
183
184 #define PAR ((struct unur_gibbs_par*)par->datap) /* data for parameter object */
185 #define GEN ((struct unur_gibbs_gen*)gen->datap) /* data for generator object */
186 #define DISTR gen->distr->data.cvec /* data for distribution in generator object */
187
188 #define SAMPLE gen->sample.cvec /* pointer to sampling routine */
189
190 /* generators for conditional distributions */
191 #define GEN_CONDI gen->gen_aux_list
192
193 /* an auxiliary generator for standard normal variates */
194 #define GEN_NORMAL gen->gen_aux
195
196 /*---------------------------------------------------------------------------*/
197
198 static UNUR_SAMPLING_ROUTINE_CVEC *
_unur_gibbs_getSAMPLE(struct unur_gen * gen)199 _unur_gibbs_getSAMPLE( struct unur_gen *gen )
200 {
201 switch (gen->variant & GIBBS_VARMASK_VARIANT) {
202 case GIBBS_VARIANT_RANDOMDIR:
203 return _unur_gibbs_randomdir_sample_cvec;
204 case GIBBS_VARIANT_COORD:
205 default:
206 return _unur_gibbs_coord_sample_cvec;
207 }
208 } /* end of _unur_gibbs_getSAMPLE() */
209
210 /*---------------------------------------------------------------------------*/
211
212 /*****************************************************************************/
213 /** Public: User Interface (API) **/
214 /*****************************************************************************/
215
216 struct unur_par *
unur_gibbs_new(const struct unur_distr * distr)217 unur_gibbs_new( const struct unur_distr *distr )
218 /*----------------------------------------------------------------------*/
219 /* get default parameters */
220 /* */
221 /* parameters: */
222 /* distr ... pointer to distribution object */
223 /* */
224 /* return: */
225 /* default parameters (pointer to structure) */
226 /* */
227 /* error: */
228 /* return NULL */
229 /*----------------------------------------------------------------------*/
230 {
231 struct unur_par *par;
232
233 /* check arguments */
234 _unur_check_NULL( GENTYPE,distr,NULL );
235
236 /* check distribution */
237 if (distr->type != UNUR_DISTR_CVEC) {
238 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
239 COOKIE_CHECK(distr,CK_DISTR_CVEC,NULL);
240
241 if (DISTR_IN.logpdf == NULL) {
242 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"logPDF");
243 return NULL;
244 }
245 if (DISTR_IN.dlogpdf == NULL) {
246 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"dlogPDF");
247 return NULL;
248 }
249
250 /* allocate structure */
251 par = _unur_par_new( sizeof(struct unur_gibbs_par) );
252 COOKIE_SET(par,CK_GIBBS_PAR);
253
254 /* copy input */
255 par->distr = distr; /* pointer to distribution object */
256
257 /* set default values */
258 PAR->c_T = 0.; /* parameter for transformation (-1. <= c < 0.) */
259
260 par->method = UNUR_METH_GIBBS ; /* method */
261 par->variant = GIBBS_VARIANT_COORD; /* default variant */
262 par->set = 0u; /* inidicate default parameters */
263 par->urng = unur_get_default_urng(); /* use default urng */
264 par->urng_aux = NULL; /* no auxilliary URNG required */
265
266 PAR->thinning = 1; /* thinning factor */
267 PAR->burnin = 0; /* length of burn-in for chain */
268 PAR->x0 = NULL; /* starting point of chain, default is 0 */
269
270 par->debug = _unur_default_debugflag; /* set default debugging flags */
271
272 /* routine for starting generator */
273 par->init = _unur_gibbs_init;
274
275 return par;
276
277 } /* end of unur_gibbs_new() */
278
279 /*****************************************************************************/
280
281 int
unur_gibbs_set_variant_coordinate(struct unur_par * par)282 unur_gibbs_set_variant_coordinate( struct unur_par *par )
283 /*----------------------------------------------------------------------*/
284 /* Coordinate Sampler : */
285 /* Sampling along the coordinate directions (cyclic). */
286 /* */
287 /* parameters: */
288 /* par ... pointer to parameter for building generator object */
289 /* */
290 /* return: */
291 /* UNUR_SUCCESS ... on success */
292 /* error code ... on error */
293 /*----------------------------------------------------------------------*/
294 {
295 /* check arguments */
296 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
297 _unur_check_par_object( par, GIBBS );
298
299 /* we use a bit in variant */
300 par->variant = (par->variant & ~GIBBS_VARMASK_VARIANT) | GIBBS_VARIANT_COORD;
301
302 /* ok */
303 return UNUR_SUCCESS;
304 } /* end of unur_gibbs_set_variant_coordinate() */
305
306 /*---------------------------------------------------------------------------*/
307
308 int
unur_gibbs_set_variant_random_direction(struct unur_par * par)309 unur_gibbs_set_variant_random_direction( struct unur_par *par )
310 /*----------------------------------------------------------------------*/
311 /* Random Direction Sampler : */
312 /* Sampling along the random directions. */
313 /* */
314 /* parameters: */
315 /* par ... pointer to parameter for building generator object */
316 /* */
317 /* return: */
318 /* UNUR_SUCCESS ... on success */
319 /* error code ... on error */
320 /*----------------------------------------------------------------------*/
321 {
322 /* check arguments */
323 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
324 _unur_check_par_object( par, GIBBS );
325
326 /* we use a bit in variant */
327 par->variant = (par->variant & ~GIBBS_VARMASK_VARIANT) | GIBBS_VARIANT_RANDOMDIR;
328
329 /* ok */
330 return UNUR_SUCCESS;
331 } /* end of unur_gibbs_set_variant_coordinate() */
332
333 /*---------------------------------------------------------------------------*/
334
335 int
unur_gibbs_set_c(struct unur_par * par,double c)336 unur_gibbs_set_c( struct unur_par *par, double c )
337 /*----------------------------------------------------------------------*/
338 /* set parameter c for transformation T_c used for sampling from */
339 /* conditional distributions */
340 /* */
341 /* parameters: */
342 /* par ... pointer to parameter for building generator object */
343 /* c ... parameter c */
344 /* */
345 /* return: */
346 /* UNUR_SUCCESS ... on success */
347 /* error code ... on error */
348 /*----------------------------------------------------------------------*/
349 {
350 /* check arguments */
351 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
352 _unur_check_par_object( par, GIBBS );
353
354 /* check new parameter for generator */
355 if (c > 0.) {
356 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"c > 0");
357 return UNUR_ERR_PAR_SET;
358 }
359 /** TODO: ... **/
360 /* if (c <= -1.) { */
361 /* _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"c <= -1 only if domain is bounded. Use `TABL' method then."); */
362 /* return 0; */
363 /* } */
364 /** TODO: ... **/
365 if (c < -0.5) {
366 _unur_error(GENTYPE,UNUR_ERR_PAR_SET,"c < -0.5 not implemented yet");
367 return UNUR_ERR_PAR_SET;
368 }
369 if (!_unur_iszero(c) && c > -0.5) {
370 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"-0.5 < c < 0 not recommended. using c = -0.5 instead.");
371 c = -0.5;
372 }
373
374 /* store date */
375 PAR->c_T = c;
376
377 /* changelog */
378 par->set |= GIBBS_SET_C;
379
380 return UNUR_SUCCESS;
381
382 } /* end of unur_gibbs_set_c() */
383
384 /*---------------------------------------------------------------------------*/
385
386 int
unur_gibbs_set_startingpoint(struct unur_par * par,const double * x0)387 unur_gibbs_set_startingpoint( struct unur_par *par, const double *x0)
388 /*----------------------------------------------------------------------*/
389 /* set starting point for chain */
390 /* */
391 /* parameters: */
392 /* par ... pointer to parameter for building generator object */
393 /* x0 ... starting point of chain */
394 /* */
395 /* return: */
396 /* UNUR_SUCCESS ... on success */
397 /* error code ... on error */
398 /*----------------------------------------------------------------------*/
399 {
400 /* check arguments */
401 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
402 _unur_check_par_object( par, GIBBS );
403
404 /* store data */
405 PAR->x0 = x0;
406
407 /* changelog */
408 par->set |= GIBBS_SET_X0;
409
410 /* ok */
411 return UNUR_SUCCESS;
412 } /* end of unur_gibbs_set_startingpoint() */
413
414 /*---------------------------------------------------------------------------*/
415
416 int
unur_gibbs_set_thinning(struct unur_par * par,int thinning)417 unur_gibbs_set_thinning( struct unur_par *par, int thinning )
418 /*----------------------------------------------------------------------*/
419 /* set thinning factor for chain */
420 /* */
421 /* parameters: */
422 /* par ... pointer to parameter for building generator object */
423 /* thinning ... thinning factor */
424 /* */
425 /* return: */
426 /* UNUR_SUCCESS ... on success */
427 /* error code ... on error */
428 /*----------------------------------------------------------------------*/
429 {
430 /* check arguments */
431 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
432 _unur_check_par_object( par, GIBBS );
433
434 /* check new parameter for generator */
435 if (thinning < 1) {
436 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"thinning < 1");
437 return UNUR_ERR_PAR_SET;
438 }
439
440 /* store data */
441 PAR->thinning = thinning;
442
443 /* changelog */
444 par->set |= GIBBS_SET_THINNING;
445
446 /* o.k. */
447 return UNUR_SUCCESS;
448
449 } /* end of unur_gibbs_set_thinning() */
450
451 /*---------------------------------------------------------------------------*/
452
453 int
unur_gibbs_set_burnin(struct unur_par * par,int burnin)454 unur_gibbs_set_burnin( struct unur_par *par, int burnin )
455 /*----------------------------------------------------------------------*/
456 /* set length of burn-in for chain */
457 /* */
458 /* parameters: */
459 /* par ... pointer to parameter for building generator object */
460 /* burnin ... length of burn-in */
461 /* */
462 /* return: */
463 /* UNUR_SUCCESS ... on success */
464 /* error code ... on error */
465 /*----------------------------------------------------------------------*/
466 {
467 /* check arguments */
468 _unur_check_NULL( GENTYPE, par, UNUR_ERR_NULL );
469 _unur_check_par_object( par, GIBBS );
470
471 /* check new parameter for generator */
472 if (burnin < 0) {
473 _unur_warning(GENTYPE,UNUR_ERR_PAR_SET,"burnin < 0");
474 return UNUR_ERR_PAR_SET;
475 }
476
477 /* store data */
478 PAR->burnin = burnin;
479
480 /* changelog */
481 par->set |= GIBBS_SET_BURNIN;
482
483 /* o.k. */
484 return UNUR_SUCCESS;
485
486 } /* end of unur_gibbs_set_burnin() */
487
488 /*---------------------------------------------------------------------------*/
489
490 const double *
unur_gibbs_get_state(struct unur_gen * gen)491 unur_gibbs_get_state( struct unur_gen *gen )
492 /*----------------------------------------------------------------------*/
493 /* get current state of the Gibbs chain */
494 /* */
495 /* parameters: */
496 /* gen ... pointer to generator object */
497 /* */
498 /* return: */
499 /* pointer to chain ... on success */
500 /* NULL ... on error */
501 /*----------------------------------------------------------------------*/
502 {
503 /* check arguments */
504 _unur_check_NULL( GENTYPE, gen, NULL );
505 if (gen->method != UNUR_METH_GIBBS) {
506 _unur_error(gen->genid, UNUR_ERR_GEN_INVALID,"");
507 return NULL;
508 }
509
510 return GEN->state;
511 } /* end of unur_gibbs_get_state() */
512
513 /*---------------------------------------------------------------------------*/
514
515 int
unur_gibbs_chg_state(struct unur_gen * gen,const double * state)516 unur_gibbs_chg_state( struct unur_gen *gen, const double *state )
517 /*----------------------------------------------------------------------*/
518 /* chg current state of the Gibbs chain */
519 /* */
520 /* parameters: */
521 /* gen ... pointer to generator object */
522 /* state ... new state of chain */
523 /* */
524 /* return: */
525 /* UNUR_SUCCESS ... on success */
526 /* error code ... on error */
527 /*----------------------------------------------------------------------*/
528 {
529 /* check arguments */
530 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
531 _unur_check_gen_object( gen, GIBBS, UNUR_ERR_GEN_INVALID );
532 _unur_check_NULL( gen->genid, state, UNUR_ERR_NULL );
533
534 /* copy state */
535 memcpy( GEN->state, state, GEN->dim * sizeof(double));
536
537 return UNUR_SUCCESS;
538 } /* end of unur_gibbs_chg_state() */
539
540 /*---------------------------------------------------------------------------*/
541
542 int
unur_gibbs_reset_state(struct unur_gen * gen)543 unur_gibbs_reset_state( struct unur_gen *gen )
544 /*----------------------------------------------------------------------*/
545 /* reset current state of the Gibbs chain to starting point */
546 /* */
547 /* parameters: */
548 /* gen ... pointer to generator object */
549 /* */
550 /* return: */
551 /* UNUR_SUCCESS ... on success */
552 /* error code ... on error */
553 /*----------------------------------------------------------------------*/
554 {
555 /* check arguments */
556 _unur_check_NULL( GENTYPE, gen, UNUR_ERR_NULL );
557 _unur_check_gen_object( gen, GIBBS, UNUR_ERR_GEN_INVALID );
558
559 /* copy state */
560 memcpy( GEN->state, GEN->x0, GEN->dim * sizeof(double));
561
562 if (gen->variant & GIBBS_VARIANT_COORD)
563 GEN->coord = (GEN->dim)-1;
564
565 return UNUR_SUCCESS;
566 } /* end of unur_gibbs_reset_state() */
567
568
569 /*****************************************************************************/
570 /** Private **/
571 /*****************************************************************************/
572
573 struct unur_gen *
_unur_gibbs_init(struct unur_par * par)574 _unur_gibbs_init( struct unur_par *par )
575 /*----------------------------------------------------------------------*/
576 /* initialize new generator */
577 /* */
578 /* parameters: */
579 /* par ... pointer to paramters for building generator object */
580 /* */
581 /* return: */
582 /* pointer to generator object */
583 /* */
584 /* error: */
585 /* return NULL */
586 /*----------------------------------------------------------------------*/
587 {
588 struct unur_gen *gen;
589
590 /* check arguments */
591 _unur_check_NULL( GENTYPE,par,NULL );
592
593 /* check input */
594 if ( par->method != UNUR_METH_GIBBS ) {
595 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
596 return NULL; }
597 COOKIE_CHECK(par,CK_GIBBS_PAR,NULL);
598
599 /* create a new empty generator object */
600 gen = _unur_gibbs_create(par);
601 _unur_par_free(par);
602 if (!gen) return NULL;
603
604 #ifdef UNUR_ENABLE_LOGGING
605 /* write info into LOG file */
606 if (gen->debug) _unur_gibbs_debug_init_start(gen);
607 #endif
608
609 /* make generators for conditional distributions */
610 switch (gen->variant & GIBBS_VARMASK_VARIANT) {
611 case GIBBS_VARIANT_COORD: /* --- coordinate direction sampling --- */
612 if (_unur_gibbs_coord_init(gen)!=UNUR_SUCCESS) {
613 _unur_gibbs_free(gen); return NULL;
614 }
615 break;
616 case GIBBS_VARIANT_RANDOMDIR: /* --- random direction sampling --- */
617 if (_unur_gibbs_randomdir_init(gen)!=UNUR_SUCCESS) {
618 _unur_gibbs_free(gen); return NULL;
619 }
620 break;
621 default:
622 _unur_error(GENTYPE,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
623 _unur_gibbs_free(gen); return NULL;
624 }
625
626 #ifdef UNUR_ENABLE_LOGGING
627 if (gen->debug) _unur_gibbs_debug_init_condi(gen);
628 #endif
629
630 /* run burn-in */
631 if (GEN->burnin > 0 ) {
632 int thinning, burnin;
633 double *X;
634
635 /* allocate memory for random vector */
636 X = _unur_xmalloc( GEN->dim * sizeof(double) );
637
638 /* store thinning factor; we use 1 for burn-in */
639 thinning = GEN->thinning;
640 GEN->thinning = 1;
641
642 for (burnin = GEN->burnin; burnin>0; --burnin) {
643 if ( _unur_sample_vec(gen,X) != UNUR_SUCCESS ) {
644 #ifdef UNUR_ENABLE_LOGGING
645 _unur_gibbs_debug_burnin_failed(gen);
646 if (gen->debug) _unur_gibbs_debug_init_finished(gen,FALSE);
647 #endif
648 _unur_gibbs_free(gen); free (X); return NULL;
649 }
650 }
651
652 /* restore thinning factor */
653 GEN->thinning = thinning;
654 /* free memory for random vector */
655 free (X);
656 }
657
658 #ifdef UNUR_ENABLE_LOGGING
659 /* write info into LOG file */
660 if (gen->debug) _unur_gibbs_debug_init_finished(gen,TRUE);
661 #endif
662
663 /* o.k. */
664 return gen;
665
666 } /* end of _unur_gibbs_init() */
667
668 /*---------------------------------------------------------------------------*/
669
670 int
_unur_gibbs_coord_init(struct unur_gen * gen)671 _unur_gibbs_coord_init( struct unur_gen *gen )
672 /*----------------------------------------------------------------------*/
673 /* initialize new generator for coordinate sampler */
674 /* */
675 /* parameters: */
676 /* gen ... pointer to generator object */
677 /* */
678 /* return: */
679 /* UNUR_SUCCESS ... on success */
680 /* error code ... on error */
681 /*----------------------------------------------------------------------*/
682 {
683 struct unur_par *par_condi;
684 struct unur_gen *gen_condi;
685 int i;
686 int errorcode = UNUR_SUCCESS;
687
688 /* conditional distribution object */
689 GEN->distr_condi = unur_distr_condi_new( gen->distr, GEN->state, NULL, 0);
690
691 for (i=0; i<GEN->dim; i++) {
692
693 if ( (errorcode=unur_distr_condi_set_condition(GEN->distr_condi,GEN->state,NULL,i))
694 != UNUR_SUCCESS )
695 break;
696
697 /* make parameter object */
698 switch( gen->variant & GIBBS_VARMASK_T ) {
699 case GIBBS_VAR_T_LOG:
700 /* use more robust method ARS for T = log */
701 par_condi = unur_ars_new(GEN->distr_condi);
702 unur_ars_set_reinit_percentiles(par_condi,2,NULL);
703 break;
704
705 case GIBBS_VAR_T_SQRT:
706 /* we only have method TDR for T = -1/sqrt */
707 par_condi = unur_tdr_new(GEN->distr_condi);
708 unur_tdr_set_reinit_percentiles(par_condi,2,NULL);
709 unur_tdr_set_c(par_condi,-0.5);
710 unur_tdr_set_usedars(par_condi,FALSE);
711 unur_tdr_set_variant_gw(par_condi);
712 break;
713
714 case GIBBS_VAR_T_POW:
715 default:
716 _unur_error(gen->genid,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
717 return UNUR_ERR_SHOULD_NOT_HAPPEN;
718 }
719
720 /* we do not need a private copy since otherwise we cannot update the generator object */
721 unur_set_use_distr_privatecopy( par_condi, FALSE );
722 /* debugging messages from the conditional generator should be rare */
723 /* otherwise the size if the LOG file explodes */
724 unur_set_debug( par_condi, (gen->debug&GIBBS_DEBUG_CONDI)?gen->debug:(gen->debug?1u:0u));
725
726 /* we use the same URNG for all auxiliary generators */
727 unur_set_urng( par_condi, gen->urng );
728
729 /* init generator object for sampling from conditional distributions */
730 gen_condi = unur_init(par_condi);
731
732 if (gen_condi == NULL) {
733 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,
734 "Cannot create generator for conditional distributions");
735 #ifdef UNUR_ENABLE_LOGGING
736 if (gen->debug) _unur_gibbs_debug_init_finished(gen,FALSE);
737 #endif
738 errorcode = UNUR_ERR_GEN_CONDITION;
739 break;
740 }
741
742 /* store pointer to generator of conditional distribution */
743 GEN_CONDI[i] = gen_condi;
744
745 if (i==0 && DISTR.domainrect==NULL) {
746 /* when the domain is not bounded by a rectangle we only */
747 /* need a clone for each dimension (except the first one) */
748 for (i=1; i<GEN->dim; i++)
749 GEN_CONDI[i] = unur_gen_clone(gen_condi);
750 break;
751 }
752 }
753
754 return errorcode;
755 } /* end of _unur_gibbs_coord_init() */
756
757 /*---------------------------------------------------------------------------*/
758
759 int
_unur_gibbs_randomdir_init(struct unur_gen * gen)760 _unur_gibbs_randomdir_init( struct unur_gen *gen )
761 /*----------------------------------------------------------------------*/
762 /* initialize new generator for random direction sampler */
763 /* */
764 /* parameters: */
765 /* gen ... pointer to generator object */
766 /* */
767 /* return: */
768 /* UNUR_SUCCESS ... on success */
769 /* error code ... on error */
770 /*----------------------------------------------------------------------*/
771 {
772 struct unur_par *par_condi;
773 struct unur_gen *gen_condi;
774
775 /* we need an auxiliary generator for normal random variates */
776 GEN_NORMAL = _unur_gibbs_normalgen( gen );
777 if ( GEN_NORMAL == NULL ) return UNUR_FAILURE;
778
779 /* conditional distribution object */
780 _unur_gibbs_random_unitvector( gen, GEN->direction );
781 GEN->distr_condi = unur_distr_condi_new( gen->distr, GEN->state, GEN->direction, 0);
782
783 /* make parameter object */
784 switch( gen->variant & GIBBS_VARMASK_T ) {
785 case GIBBS_VAR_T_LOG:
786 /* use more robust method ARS for T = log */
787 par_condi = unur_ars_new(GEN->distr_condi);
788 unur_ars_set_reinit_percentiles(par_condi,2,NULL);
789 break;
790
791 case GIBBS_VAR_T_SQRT:
792 /* we only have method TDR for T = -1/sqrt */
793 par_condi = unur_tdr_new(GEN->distr_condi);
794 unur_tdr_set_reinit_percentiles(par_condi,2,NULL);
795 unur_tdr_set_c(par_condi,-0.5);
796 unur_tdr_set_usedars(par_condi,FALSE);
797 break;
798
799 case GIBBS_VAR_T_POW:
800 default:
801 _unur_error(gen->genid,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
802 return UNUR_ERR_SHOULD_NOT_HAPPEN;
803 }
804
805 /* we do not need a private copy since otherwise we cannot update the generator object */
806 unur_set_use_distr_privatecopy( par_condi, FALSE );
807 /* debugging messages from the conditional generator should be rare */
808 /* otherwise the size if the LOG file explodes */
809 unur_set_debug( par_condi, (gen->debug&GIBBS_DEBUG_CONDI)?gen->debug:(gen->debug?1u:0u));
810
811 /* we use the same URNG for all auxiliary generators */
812 unur_set_urng( par_condi, gen->urng );
813
814 /* init generator object for sampling from conditional distributions */
815 gen_condi = unur_init(par_condi);
816
817 if (gen_condi == NULL) {
818 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,
819 "Cannot create generator for conditional distributions");
820 #ifdef UNUR_ENABLE_LOGGING
821 if (gen->debug) _unur_gibbs_debug_init_finished(gen,FALSE);
822 #endif
823 return UNUR_ERR_GEN_CONDITION;
824 }
825
826 /* store generator in structure. we only need one such generator */
827 *GEN_CONDI = gen_condi;
828
829 return UNUR_SUCCESS;
830 } /* end of _unur_gibbs_randomdir_init() */
831
832 /*---------------------------------------------------------------------------*/
833
834 struct unur_gen *
_unur_gibbs_create(struct unur_par * par)835 _unur_gibbs_create( struct unur_par *par )
836 /*----------------------------------------------------------------------*/
837 /* allocate memory for generator */
838 /* */
839 /* parameters: */
840 /* par ... pointer to parameter for building generator object */
841 /* */
842 /* return: */
843 /* pointer to (empty) generator object with default settings */
844 /* */
845 /* error: */
846 /* return NULL */
847 /*----------------------------------------------------------------------*/
848 {
849 struct unur_gen *gen;
850 int i;
851
852 /* check arguments */
853 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_GIBBS_PAR,NULL);
854
855 /* create new generic generator object */
856 gen = _unur_generic_create( par, sizeof(struct unur_gibbs_gen) );
857
858 /* magic cookies */
859 COOKIE_SET(gen,CK_GIBBS_GEN);
860
861 /* dimension of distribution */
862 GEN->dim = gen->distr->dim;
863
864 /* set generator identifier */
865 gen->genid = _unur_set_genid(GENTYPE);
866
867 /* which transformation for conditional distributions */
868 if ( _unur_iszero(PAR->c_T) )
869 par->variant = (par->variant & (~GIBBS_VARMASK_T)) | GIBBS_VAR_T_LOG;
870 else if (_unur_FP_same(PAR->c_T, -0.5))
871 par->variant = (par->variant & (~GIBBS_VARMASK_T)) | GIBBS_VAR_T_SQRT;
872 else
873 par->variant = (par->variant & (~GIBBS_VARMASK_T)) | GIBBS_VAR_T_POW;
874
875 /* routines for sampling and destroying generator */
876 SAMPLE = _unur_gibbs_getSAMPLE(gen);
877 gen->destroy = _unur_gibbs_free;
878 gen->clone = _unur_gibbs_clone;
879
880 /* variant of sampling method */
881 gen->variant = par->variant;
882
883 /* copy parameters into generator object */
884 GEN->thinning = PAR->thinning; /* thinning factor */
885 GEN->burnin = PAR->burnin; /* length of burnin */
886 GEN->c_T = PAR->c_T; /* parameter for transformation */
887
888 /* allocate memory for state */
889 GEN->state = _unur_xmalloc( GEN->dim * sizeof(double));
890 GEN->x0 = _unur_xmalloc( GEN->dim * sizeof(double));
891 if (PAR->x0 == NULL)
892 PAR->x0 = unur_distr_cvec_get_center(gen->distr);
893 memcpy( GEN->state, PAR->x0, GEN->dim * sizeof(double));
894 memcpy( GEN->x0, PAR->x0, GEN->dim * sizeof(double));
895
896 /* generator(s) for conditional distributions */
897 GEN->distr_condi = NULL;
898 GEN_CONDI = _unur_xmalloc( GEN->dim * sizeof(struct unur_gen *) );
899 gen->n_gen_aux_list = GEN->dim; /* size of array GEN_CONDI */
900 for (i=0; i<GEN->dim; i++) GEN_CONDI[i] = NULL;
901
902 /* allocate memory for random direction */
903 GEN->direction = _unur_xmalloc( GEN->dim * sizeof(double));
904
905 /* defaults */
906 GEN->coord = (GEN->dim)-1; /* current coordinate of GIBBS chain.
907 we want to start with coordinate 0. */
908
909 #ifdef UNUR_ENABLE_INFO
910 /* set function for creating info string */
911 gen->info = _unur_gibbs_info;
912 #endif
913
914 /* return pointer to (almost empty) generator object */
915 return gen;
916
917 } /* end of _unur_gibbs_create() */
918
919 /*---------------------------------------------------------------------------*/
920
921 struct unur_gen *
_unur_gibbs_clone(const struct unur_gen * gen)922 _unur_gibbs_clone( const struct unur_gen *gen )
923 /*----------------------------------------------------------------------*/
924 /* copy (clone) generator object */
925 /* */
926 /* parameters: */
927 /* gen ... pointer to generator object */
928 /* */
929 /* return: */
930 /* pointer to clone of generator object */
931 /* */
932 /* error: */
933 /* return NULL */
934 /*----------------------------------------------------------------------*/
935 {
936 #define CLONE ((struct unur_gibbs_gen*)clone->datap)
937 #define CLONE_CONDI clone->gen_aux_list
938
939 int i;
940 struct unur_gen *clone;
941
942 /* check arguments */
943 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_GIBBS_GEN,NULL);
944
945 /* create generic clone */
946 clone = _unur_generic_clone( gen, GENTYPE );
947
948 /* copy state */
949 CLONE->state = _unur_xmalloc( GEN->dim * sizeof(double));
950 memcpy( CLONE->state, GEN->state, GEN->dim * sizeof(double));
951 CLONE->x0 = _unur_xmalloc( GEN->dim * sizeof(double));
952 memcpy( CLONE->x0, GEN->x0, GEN->dim * sizeof(double));
953
954 /* copy generators for conditional distributions */
955 if (GEN->distr_condi) CLONE->distr_condi = _unur_distr_clone( GEN->distr_condi );
956
957 /* GEN_CONDI is cloned by _unur_generic_clone */
958 /* however, these use a pointer to GEN->distr which must be updated, too */
959 if (CLONE_CONDI) {
960 for (i=0; i<GEN->dim; i++)
961 if (CLONE_CONDI[i])
962 CLONE_CONDI[i]->distr = CLONE->distr_condi;
963 }
964
965 /* allocate memory for random direction */
966 CLONE->direction = _unur_xmalloc( GEN->dim * sizeof(double));
967
968 return clone;
969
970 #undef CLONE
971 #undef CLONE_CONDI
972 } /* end of _unur_gibbs_clone() */
973
974 /*---------------------------------------------------------------------------*/
975
976 void
_unur_gibbs_free(struct unur_gen * gen)977 _unur_gibbs_free( struct unur_gen *gen )
978 /*----------------------------------------------------------------------*/
979 /* deallocate generator object */
980 /* */
981 /* parameters: */
982 /* gen ... pointer to generator object */
983 /*----------------------------------------------------------------------*/
984 {
985 /* check arguments */
986 if( !gen ) /* nothing to do */
987 return;
988
989 /* check input */
990 if ( gen->method != UNUR_METH_GIBBS ) {
991 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
992 return; }
993 COOKIE_CHECK(gen,CK_GIBBS_GEN,RETURN_VOID);
994
995 /* we cannot use this generator object any more */
996 SAMPLE = NULL; /* make sure to show up a programming error */
997
998 /* free vectors */
999 if (GEN->state) free (GEN->state);
1000 if (GEN->x0) free (GEN->x0);
1001 if (GEN->direction) free (GEN->direction);
1002
1003 /* free conditional distribution object */
1004 if (GEN->distr_condi) _unur_distr_free (GEN->distr_condi);
1005
1006 /* GEN_CONDI is freed by _unur_generic_free */
1007
1008 _unur_generic_free(gen);
1009
1010 } /* end of _unur_gibbs_free() */
1011
1012 /*****************************************************************************/
1013
1014 int
_unur_gibbs_coord_sample_cvec(struct unur_gen * gen,double * vec)1015 _unur_gibbs_coord_sample_cvec( struct unur_gen *gen, double *vec )
1016 /*----------------------------------------------------------------------*/
1017 /* sample from generator */
1018 /* */
1019 /* parameters: */
1020 /* gen ... pointer to generator object */
1021 /* vec ... random vector (result) */
1022 /*----------------------------------------------------------------------*/
1023 {
1024 double X;
1025 int thinning;
1026
1027 /* check arguments */
1028 CHECK_NULL(gen,UNUR_ERR_NULL);
1029 COOKIE_CHECK(gen,CK_GIBBS_GEN,UNUR_ERR_COOKIE);
1030
1031 for (thinning = GEN->thinning; thinning > 0; --thinning) {
1032
1033 /* update coordinate direction */
1034 GEN->coord = (GEN->coord + 1) % GEN->dim;
1035
1036 /* check state of chain */
1037 if (!_unur_isfinite(GEN->state[GEN->coord]))
1038 /* there has been a fatal error during the last sampling from
1039 the conditional distribution. This probably has been caused
1040 by improper target distribution for which the method
1041 does not work (i.e., it is not T-concave).
1042 */
1043 continue;
1044
1045 /* update conditional distribution */
1046 unur_distr_condi_set_condition( GEN->distr_condi, GEN->state, NULL, GEN->coord);
1047
1048 /* reinit generator object */
1049 if (unur_reinit(GEN_CONDI[GEN->coord]) == UNUR_SUCCESS) {
1050
1051 /* sample from distribution */
1052 X = unur_sample_cont(GEN_CONDI[GEN->coord]);
1053
1054 if (_unur_isfinite(X)) {
1055 /* update state */
1056 GEN->state[GEN->coord] = X;
1057 continue;
1058 }
1059 }
1060
1061 /* ERROR: reset to starting point */
1062 _unur_warning(gen->genid,UNUR_ERR_GEN_SAMPLING,"reset chain");
1063 unur_gibbs_reset_state(gen);
1064 return UNUR_FAILURE;
1065
1066 }
1067
1068 /* copy current state into given vector */
1069 memcpy(vec, GEN->state, GEN->dim * sizeof(double));
1070
1071 return UNUR_SUCCESS;
1072
1073 } /* end of _unur_gibbs_coord_sample_cvec() */
1074
1075 /*---------------------------------------------------------------------------*/
1076
1077 int
_unur_gibbs_randomdir_sample_cvec(struct unur_gen * gen,double * vec)1078 _unur_gibbs_randomdir_sample_cvec( struct unur_gen *gen, double *vec )
1079 /*----------------------------------------------------------------------*/
1080 /* sample from generator */
1081 /* */
1082 /* parameters: */
1083 /* gen ... pointer to generator object */
1084 /* vec ... random vector (result) */
1085 /*----------------------------------------------------------------------*/
1086 {
1087 int i;
1088 double X;
1089 int thinning;
1090
1091 /* check arguments */
1092 CHECK_NULL(gen,UNUR_ERR_NULL);
1093 COOKIE_CHECK(gen,CK_GIBBS_GEN,UNUR_ERR_COOKIE);
1094
1095 for (thinning = GEN->thinning; thinning > 0; --thinning) {
1096
1097 /* check state of chain */
1098 if (!_unur_isfinite(GEN->state[0]))
1099 /* there has been a fatal error during the last sampling from
1100 the conditional distribution. This probably has been caused
1101 by improper target distribution for which the method
1102 does not work (i.e., it is not T-concave).
1103 */
1104 break;
1105
1106 /* new random direction */
1107 _unur_gibbs_random_unitvector( gen, GEN->direction );
1108
1109 /* update conditional distribution */
1110 unur_distr_condi_set_condition( GEN->distr_condi, GEN->state, GEN->direction, 0);
1111
1112 /* reinit generator object */
1113 if (unur_reinit(*GEN_CONDI) == UNUR_SUCCESS) {
1114
1115 /* sample from distribution */
1116 X = unur_sample_cont(*GEN_CONDI);
1117
1118 if (_unur_isfinite(X)) {
1119 /* update state */
1120 for (i=0; i<GEN->dim; i++)
1121 GEN->state[i] += X * GEN->direction[i];
1122 continue;
1123 }
1124 }
1125
1126 /* ERROR: reset to starting point */
1127 _unur_warning(gen->genid,UNUR_ERR_GEN_SAMPLING,"reset chain");
1128 unur_gibbs_reset_state(gen);
1129 return UNUR_FAILURE;
1130
1131 }
1132
1133 /* copy current state into given vector */
1134 memcpy(vec, GEN->state, GEN->dim * sizeof(double));
1135
1136 return UNUR_SUCCESS;
1137
1138 } /* end of _unur_gibbs_randomdir_sample_cvec() */
1139
1140
1141 /*****************************************************************************/
1142 /** Auxilliary Routines **/
1143 /*****************************************************************************/
1144
1145 struct unur_gen *
_unur_gibbs_normalgen(struct unur_gen * gen)1146 _unur_gibbs_normalgen( struct unur_gen *gen )
1147 /*----------------------------------------------------------------------*/
1148 /* create a normal random variate generator */
1149 /* */
1150 /* parameters: */
1151 /* gen ... pointer to GIBBS generator object */
1152 /* */
1153 /* */
1154 /*----------------------------------------------------------------------*/
1155 {
1156 struct unur_gen *normalgen;
1157 struct unur_distr *normaldistr = unur_distr_normal(NULL,0);
1158 struct unur_par *normalpar = unur_arou_new( normaldistr );
1159
1160 unur_arou_set_usedars( normalpar, TRUE );
1161 normalgen = unur_init( normalpar );
1162 _unur_distr_free( normaldistr );
1163 if (normalgen == NULL) {
1164 _unur_error(gen->genid,UNUR_ERR_SHOULD_NOT_HAPPEN,
1165 "Cannot create aux Gaussian generator");
1166 return NULL;
1167 }
1168
1169 /* uniform random number generator and debugging flags */
1170 normalgen->urng = gen->urng;
1171 normalgen->debug = gen->debug;
1172
1173 return normalgen;
1174
1175 } /* end of _unur_gibbs_normalgen() */
1176
1177 /*---------------------------------------------------------------------------*/
1178
1179 void
_unur_gibbs_random_unitvector(struct unur_gen * gen,double * direction)1180 _unur_gibbs_random_unitvector( struct unur_gen *gen, double *direction )
1181 /*----------------------------------------------------------------------*/
1182 /* generate a random direction vector */
1183 /* */
1184 /* parameters: */
1185 /* gen ... pointer to generator object */
1186 /* direction ... random vector (result) */
1187 /*----------------------------------------------------------------------*/
1188 {
1189 int i;
1190
1191 do {
1192 for (i=0; i<GEN->dim; i++)
1193 direction[i] = unur_sample_cont(GEN_NORMAL);
1194 /* normalize direction vector */
1195 _unur_vector_normalize(GEN->dim, direction);
1196
1197 /* there is an extremely small change that direction is the null before
1198 normalizing. In this case non of its coordinates are finite. */
1199 } while (!_unur_isfinite(direction[0]));
1200
1201 } /* end of _unur_gibbs_random_unitvector() */
1202
1203
1204 /*****************************************************************************/
1205 /** Debugging utilities **/
1206 /*****************************************************************************/
1207
1208 /*---------------------------------------------------------------------------*/
1209 #ifdef UNUR_ENABLE_LOGGING
1210 /*---------------------------------------------------------------------------*/
1211
1212 void
_unur_gibbs_debug_init_start(const struct unur_gen * gen)1213 _unur_gibbs_debug_init_start( const struct unur_gen *gen )
1214 /*----------------------------------------------------------------------*/
1215 /* write info about generator into LOG file */
1216 /* */
1217 /* parameters: */
1218 /* gen ... pointer to generator object */
1219 /*----------------------------------------------------------------------*/
1220 {
1221 FILE *LOG;
1222
1223 /* check arguments */
1224 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_GIBBS_GEN,RETURN_VOID);
1225
1226 LOG = unur_get_stream();
1227
1228 fprintf(LOG,"%s:\n",gen->genid);
1229 fprintf(LOG,"%s: type = continuous multivariate random variates\n",gen->genid);
1230 fprintf(LOG,"%s: method = GIBBS (Markov Chain - GIBBS sampler)\n",gen->genid);
1231 fprintf(LOG,"%s: variant = ",gen->genid);
1232 switch (gen->variant & GIBBS_VARMASK_VARIANT) {
1233 case GIBBS_VARIANT_COORD:
1234 fprintf(LOG,"coordinate sampling (original Gibbs sampler) [default]\n"); break;
1235 case GIBBS_VARIANT_RANDOMDIR:
1236 fprintf(LOG,"random directions\n"); break;
1237 }
1238 fprintf(LOG,"%s:\n",gen->genid);
1239
1240 fprintf(LOG,"%s: transformation T_c(x) for TDR method = ",gen->genid);
1241 switch( gen->variant & GIBBS_VARMASK_T ) {
1242 case GIBBS_VAR_T_LOG:
1243 fprintf(LOG,"log(x) ... c = 0"); break;
1244 case GIBBS_VAR_T_SQRT:
1245 fprintf(LOG,"-1/sqrt(x) ... c = -1/2"); break;
1246 case GIBBS_VAR_T_POW:
1247 fprintf(LOG,"-x^(%g) ... c = %g",GEN->c_T,GEN->c_T); break;
1248 }
1249 _unur_print_if_default(gen,GIBBS_SET_C);
1250 fprintf(LOG,"\n%s:\n",gen->genid);
1251
1252 _unur_distr_cvec_debug( gen->distr, gen->genid );
1253
1254 switch (gen->variant & GIBBS_VARMASK_VARIANT) {
1255 case GIBBS_VARIANT_COORD:
1256 fprintf(LOG,"%s: sampling routine = _unur_gibbs_coord_sample()\n",gen->genid);
1257 break;
1258 case GIBBS_VARIANT_RANDOMDIR:
1259 fprintf(LOG,"%s: sampling routine = _unur_gibbs_randomdir_sample()\n",gen->genid);
1260 break;
1261 }
1262
1263 fprintf(LOG,"%s: thinning = %d",gen->genid,GEN->thinning);
1264 _unur_print_if_default(gen,GIBBS_SET_THINNING);
1265 fprintf(LOG,"\n%s: burn-in = %d",gen->genid,GEN->burnin);
1266 _unur_print_if_default(gen,GIBBS_SET_BURNIN);
1267 fprintf(LOG,"\n%s:\n",gen->genid);
1268 _unur_matrix_print_vector( GEN->dim, GEN->x0, "starting point = ", LOG, gen->genid, "\t ");
1269
1270 fprintf(LOG,"%s:\n",gen->genid);
1271
1272 } /* end of _unur_gibbs_debug_init_start() */
1273
1274 /*---------------------------------------------------------------------------*/
1275
1276 void
_unur_gibbs_debug_init_finished(const struct unur_gen * gen,int success)1277 _unur_gibbs_debug_init_finished( const struct unur_gen *gen, int success )
1278 /*----------------------------------------------------------------------*/
1279 /* write info about generator into LOG file */
1280 /* */
1281 /* parameters: */
1282 /* gen ... pointer to generator object */
1283 /* success ... whether init has failed or not */
1284 /*----------------------------------------------------------------------*/
1285 {
1286 FILE *LOG;
1287
1288 /* check arguments */
1289 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_GIBBS_GEN,RETURN_VOID);
1290
1291 LOG = unur_get_stream();
1292
1293 if (success)
1294 fprintf(LOG,"%s: INIT completed **********************\n",gen->genid);
1295 else
1296 fprintf(LOG,"%s: INIT failed **********************\n",gen->genid);
1297
1298 } /* end of _unur_gibbs_debug_init_finished() */
1299
1300 /*---------------------------------------------------------------------------*/
1301
1302 void
_unur_gibbs_debug_burnin_failed(const struct unur_gen * gen)1303 _unur_gibbs_debug_burnin_failed( const struct unur_gen *gen )
1304 /*----------------------------------------------------------------------*/
1305 /* write info after burnin has failed */
1306 /* */
1307 /* parameters: */
1308 /* gen ... pointer to generator object */
1309 /*----------------------------------------------------------------------*/
1310 {
1311 FILE *LOG;
1312
1313 /* check arguments */
1314 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_GIBBS_GEN,RETURN_VOID);
1315
1316 LOG = unur_get_stream();
1317
1318 fprintf(LOG,"%s: Burn-in failed --> INIT failed **********************\n",gen->genid);
1319 fprintf(LOG,"%s:\n",gen->genid);
1320
1321 } /* end of _unur_gibbs_debug_burnin_failed() */
1322
1323 /*---------------------------------------------------------------------------*/
1324
1325 void
_unur_gibbs_debug_init_condi(const struct unur_gen * gen)1326 _unur_gibbs_debug_init_condi( const struct unur_gen *gen )
1327 /*----------------------------------------------------------------------*/
1328 /* write list of conditional generators into LOG file */
1329 /* */
1330 /* parameters: */
1331 /* gen ... pointer to generator object */
1332 /*----------------------------------------------------------------------*/
1333 {
1334 int i;
1335 FILE *LOG;
1336
1337 /* check arguments */
1338 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_GIBBS_GEN,RETURN_VOID);
1339
1340 LOG = unur_get_stream();
1341
1342 switch (gen->variant & GIBBS_VARMASK_VARIANT) {
1343 case GIBBS_VARIANT_COORD:
1344 fprintf(LOG,"%s: generators for full conditional distributions = \n",gen->genid);
1345 fprintf(LOG,"%s:\t",gen->genid);
1346 for (i=0; i<GEN->dim; i++)
1347 fprintf(LOG,"[%s] ", GEN_CONDI[i]->genid);
1348 fprintf(LOG,"\n%s:\n",gen->genid);
1349 break;
1350 case GIBBS_VARIANT_RANDOMDIR:
1351 fprintf(LOG,"%s: generators for full conditional distributions = [%s]\n",gen->genid,
1352 GEN_CONDI[0]->genid);
1353 fprintf(LOG,"%s: generator for random directions = [%s]\n",gen->genid,
1354 GEN_NORMAL->genid);
1355 fprintf(LOG,"%s:\n",gen->genid);
1356 break;
1357 }
1358
1359 } /* end of _unur_gibbs_debug_init_condi() */
1360
1361 /*---------------------------------------------------------------------------*/
1362 #endif /* end UNUR_ENABLE_LOGGING */
1363 /*---------------------------------------------------------------------------*/
1364
1365
1366 /*---------------------------------------------------------------------------*/
1367 #ifdef UNUR_ENABLE_INFO
1368 /*---------------------------------------------------------------------------*/
1369
1370 void
_unur_gibbs_info(struct unur_gen * gen,int help)1371 _unur_gibbs_info( struct unur_gen *gen, int help )
1372 /*----------------------------------------------------------------------*/
1373 /* create character string that contains information about the */
1374 /* given generator object. */
1375 /* */
1376 /* parameters: */
1377 /* gen ... pointer to generator object */
1378 /* help ... whether to print additional comments */
1379 /*----------------------------------------------------------------------*/
1380 {
1381 struct unur_string *info = gen->infostr;
1382 struct unur_distr *distr = gen->distr;
1383 int samplesize = 10000;
1384
1385 /* generator ID */
1386 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
1387
1388 /* distribution */
1389 _unur_string_append(info,"distribution:\n");
1390 _unur_distr_info_typename(gen);
1391 _unur_string_append(info," dimension = %d\n",GEN->dim);
1392 _unur_string_append(info," functions = PDF dPDF\n");
1393 _unur_distr_cvec_info_domain(gen);
1394 _unur_string_append(info," center = ");
1395 _unur_distr_info_vector( gen, unur_distr_cvec_get_center(gen->distr), GEN->dim);
1396 if ( !(distr->set & UNUR_DISTR_SET_CENTER) ) {
1397 if ( distr->set & UNUR_DISTR_SET_MODE )
1398 _unur_string_append(info," [= mode]");
1399 else
1400 _unur_string_append(info," [default]");
1401 }
1402 _unur_string_append(info,"\n\n");
1403
1404 /* if (help) { */
1405 /* _unur_string_append(info,"\n"); */
1406 /* } */
1407
1408 /* method */
1409 _unur_string_append(info,"method: GIBBS (GIBBS sampler [MCMC])\n");
1410 _unur_string_append(info," variant = %s\n",
1411 ((gen->variant & GIBBS_VARMASK_VARIANT)==GIBBS_VARIANT_COORD)
1412 ? "coordinate sampling [default]" : "random direction sampling");
1413
1414 /* used transformation */
1415 _unur_string_append(info," T_c(x) = ");
1416 switch( gen->variant & GIBBS_VARMASK_T ) {
1417 case GIBBS_VAR_T_LOG:
1418 _unur_string_append(info,"log(x) ... c = 0\n"); break;
1419 case GIBBS_VAR_T_SQRT:
1420 _unur_string_append(info,"-1/sqrt(x) ... c = -1/2\n"); break;
1421 case GIBBS_VAR_T_POW:
1422 _unur_string_append(info,"-x^(%g) ... c = %g\n",GEN->c_T,GEN->c_T); break;
1423 }
1424
1425 _unur_string_append(info," thinning = %d\n", GEN->thinning);
1426 _unur_string_append(info,"\n");
1427
1428 /* performance */
1429 _unur_string_append(info,"performance characteristics:\n");
1430 _unur_string_append(info," rejection constant = %.2f [approx.]\n",
1431 unur_test_count_urn(gen,samplesize,0,NULL)/(2.*samplesize));
1432 _unur_string_append(info,"\n");
1433
1434 /* parameters */
1435 if (help) {
1436 _unur_string_append(info,"parameters:\n");
1437
1438 switch (gen->variant & GIBBS_VARMASK_VARIANT) {
1439 case GIBBS_VARIANT_COORD:
1440 _unur_string_append(info," variant_coordinate [default]\n"); break;
1441 case GIBBS_VARIANT_RANDOMDIR:
1442 _unur_string_append(info," variant_random_direction\n"); break;
1443 }
1444
1445 _unur_string_append(info," c = %g %s\n", GEN->c_T,
1446 (gen->set & GIBBS_SET_C) ? "" : "[default]");
1447 _unur_string_append(info," thinning = %d %s\n", GEN->thinning,
1448 (gen->set & GIBBS_SET_THINNING) ? "" : "[default]");
1449 _unur_string_append(info," burnin = %d %s\n", GEN->burnin,
1450 (gen->set & GIBBS_SET_THINNING) ? "" : "[default]");
1451
1452 _unur_string_append(info,"\n");
1453
1454 /* Not displayed:
1455 int unur_gibbs_set_startingpoint( UNUR_PAR *parameters, const double *x0 );
1456 */
1457 }
1458
1459 /* Hints */
1460 /* if (help) { */
1461 /* _unur_string_append(info,"\n"); */
1462 /* } */
1463
1464 } /* end of _unur_gibbs_info() */
1465
1466 /*---------------------------------------------------------------------------*/
1467 #endif /* end UNUR_ENABLE_INFO */
1468 /*---------------------------------------------------------------------------*/
1469