1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: norta.c *
8 * *
9 * TYPE: continuous multivariate random variate *
10 * METHOD: generates random vector with given covariance matrix and *
11 * given marginal distribution. *
12 * *
13 * DESCRIPTION: *
14 * *
15 *****************************************************************************
16 * *
17 * Copyright (c) 2000-2006 Wolfgang Hoermann and Josef Leydold *
18 * Department of Statistics and Mathematics, WU Wien, Austria *
19 * *
20 * This program is free software; you can redistribute it and/or modify *
21 * it under the terms of the GNU General Public License as published by *
22 * the Free Software Foundation; either version 2 of the License, or *
23 * (at your option) any later version. *
24 * *
25 * This program is distributed in the hope that it will be useful, *
26 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
27 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
28 * GNU General Public License for more details. *
29 * *
30 * You should have received a copy of the GNU General Public License *
31 * along with this program; if not, write to the *
32 * Free Software Foundation, Inc., *
33 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
34 * *
35 *****************************************************************************
36 * *
37 * REFERENCES: *
38 * [1] Hoermann, W., J. Leydold, and G. Derflinger (2004): *
39 * Automatic Nonuniform Random Variate Generation, Springer, Berlin. *
40 * *
41 * [2] Gosh, S. (????): Eigenvector Correction method *
42 * *
43 *****************************************************************************
44 * *
45 * NORTA (NORmal to anything) is a model to get random vectors with *
46 * given marginal distributions and rank correlation. (Notice that this *
47 * does not uniquely define the multivariate distribution.) *
48 * *
49 * In the NORTA model multinormal random variates with the given *
50 * rank (Spearman's) correlations are generated. For this task the rank *
51 * correlations are transformed into the corresponding *
52 * (Pearson) correlation matrix and the VMT method is used to sample from *
53 * the resulting multinormal distribution (by means of the Cholesky *
54 * decomposition of the covariance matrix). *
55 * In a second step the (standard normal distributed) marginal variates *
56 * are transformed by means of the CDF of the normal distribution to get *
57 * uniform marginals. The resulting random vectors then have uniform *
58 * marginals and the desired rank correlation between its components. *
59 * Such a random vector is called 'copula'. *
60 * *
61 * By means of the inverse CDF the uniform marginals are transformed into *
62 * the target marginal distributions. This transformation does not change *
63 * the rank correlation. *
64 * *
65 * It can happen that the desired rank correlation matrix is not feasible, *
66 * i.e., it cannot occur as rank correlation matrix of a multinormal *
67 * distribution. For this case we use Gosh's "eigenvector correction *
68 * method" [2]. In this case a spectral decomposition of the corresponding *
69 * (Pearson) correlation matrix is carried out. If there are some *
70 * non-positive eigenvalues (i.e. this matrix is not a true correlation *
71 * matrix) then these are set to a small value and a new (positive *
72 * definite) correlation matrix is created and used. *
73 * *
74 *****************************************************************************/
75
76 /*---------------------------------------------------------------------------*/
77
78 #include <unur_source.h>
79 #include <distr/distr.h>
80 #include <distr/distr_source.h>
81 #include <distr/cvec.h>
82 #include <distr/cont.h>
83 #include <distributions/unur_distributions.h>
84 #include <urng/urng.h>
85 #include <utils/matrix_source.h>
86 #include "unur_methods_source.h"
87 #include "x_gen.h"
88 #include "x_gen_source.h"
89 #include "cstd.h"
90 #include "hinv.h"
91 #include "ninv.h"
92 #include "pinv.h"
93 #include "mvstd.h"
94 #include "norta.h"
95 #include "norta_struct.h"
96
97 /*---------------------------------------------------------------------------*/
98 /* Constants */
99
100 /* smallest eigenvalue allowed for correlation matrix */
101 #define UNUR_NORTA_MIN_EIGENVALUE (1.e-10)
102
103 /*---------------------------------------------------------------------------*/
104 /* Variants */
105
106 /*---------------------------------------------------------------------------*/
107 /* Debugging flags */
108 /* bit 01 ... pameters and structure of generator (do not use here) */
109 /* bits 02-12 ... setup */
110 /* bits 13-24 ... adaptive steps */
111 /* bits 25-32 ... trace sampling */
112
113 #define NORTA_DEBUG_SIGMA_Y 0x00000010u /* print sigma_y for normal */
114
115 /*---------------------------------------------------------------------------*/
116 /* Flags for logging set calls */
117
118 /*---------------------------------------------------------------------------*/
119
120 #define GENTYPE "NORTA" /* type of generator */
121
122 /*---------------------------------------------------------------------------*/
123
124 static struct unur_gen *_unur_norta_init( struct unur_par *par );
125 /*---------------------------------------------------------------------------*/
126 /* Initialize new generator. */
127 /*---------------------------------------------------------------------------*/
128
129 static struct unur_gen *_unur_norta_create( struct unur_par *par );
130 /*---------------------------------------------------------------------------*/
131 /* create new (almost empty) generator object. */
132 /*---------------------------------------------------------------------------*/
133
134 static struct unur_gen *_unur_norta_clone( const struct unur_gen *gen );
135 /*---------------------------------------------------------------------------*/
136 /* copy (clone) generator object. */
137 /*---------------------------------------------------------------------------*/
138
139 static void _unur_norta_free( struct unur_gen *gen);
140 /*---------------------------------------------------------------------------*/
141 /* destroy generator object. */
142 /*---------------------------------------------------------------------------*/
143
144 static int _unur_norta_sample_cvec( struct unur_gen *gen, double *vec );
145 /*---------------------------------------------------------------------------*/
146 /* sample from generator */
147 /*---------------------------------------------------------------------------*/
148
149 static int _unur_norta_nortu_setup( struct unur_gen *gen );
150 /*---------------------------------------------------------------------------*/
151 /* Compute parameter for NORTU (normal to uniform). */
152 /*---------------------------------------------------------------------------*/
153
154 static int _unur_norta_make_correlationmatrix( int dim, double *M);
155 /*---------------------------------------------------------------------------*/
156 /* make correlation matrix by transforming symmetric positive definit matrix */
157 /*---------------------------------------------------------------------------*/
158
159 static struct unur_gen *_unur_norta_make_marginalgen( const struct unur_gen *gen,
160 const struct unur_distr *marginal );
161 /*---------------------------------------------------------------------------*/
162 /* make generator for marginal distribution */
163 /*---------------------------------------------------------------------------*/
164
165 #ifdef UNUR_ENABLE_LOGGING
166 /*---------------------------------------------------------------------------*/
167 /* the following functions print debugging information on output stream, */
168 /* i.e., into the LOG file if not specified otherwise. */
169 /*---------------------------------------------------------------------------*/
170
171 static void _unur_norta_debug_init( const struct unur_gen *gen );
172 /*---------------------------------------------------------------------------*/
173 /* print after generator has been initialized has completed. */
174 /*---------------------------------------------------------------------------*/
175
176 static void _unur_norta_debug_sigma_y( const struct unur_gen *gen,
177 const double *sigma_y,
178 const char *comment );
179 /*---------------------------------------------------------------------------*/
180 /* print sigma_y of corresponding normal distribution. */
181 /*---------------------------------------------------------------------------*/
182
183 static void _unur_norta_debug_eigensystem( const struct unur_gen *gen,
184 const double *eigenvalues,
185 const double *eigenvectors );
186 /*---------------------------------------------------------------------------*/
187 /* print eigensystem of sigma_y. */
188 /*---------------------------------------------------------------------------*/
189
190 static void _unur_norta_debug_nmgenerator( const struct unur_gen *gen );
191 /*---------------------------------------------------------------------------*/
192 /* print genid of multinormal generator. */
193 /*---------------------------------------------------------------------------*/
194 #endif
195
196 #ifdef UNUR_ENABLE_INFO
197 static void _unur_norta_info( struct unur_gen *gen, int help );
198 /*---------------------------------------------------------------------------*/
199 /* create info string. */
200 /*---------------------------------------------------------------------------*/
201 #endif
202
203 /*---------------------------------------------------------------------------*/
204 /* abbreviations */
205
206 #define DISTR_IN distr->data.cvec /* data for distribution object */
207
208 #define PAR ((struct unur_norta_par*)par->datap) /* data for parameter object */
209 #define GEN ((struct unur_norta_gen*)gen->datap) /* data for generator object */
210 #define DISTR gen->distr->data.cvec /* data for distribution in generator object */
211
212 #define SAMPLE gen->sample.cvec /* pointer to sampling routine */
213
214 #define MNORMAL gen->gen_aux /* pointer to multinormal generator */
215
216 /*---------------------------------------------------------------------------*/
217
218 #define _unur_norta_getSAMPLE(gen) (_unur_norta_sample_cvec)
219
220 /*---------------------------------------------------------------------------*/
221
222 /*****************************************************************************/
223 /** Public: User Interface (API) **/
224 /*****************************************************************************/
225
226 struct unur_par *
unur_norta_new(const struct unur_distr * distr)227 unur_norta_new( const struct unur_distr *distr )
228 /*----------------------------------------------------------------------*/
229 /* get default parameters */
230 /* */
231 /* parameters: */
232 /* distr ... pointer to distribution object */
233 /* */
234 /* return: */
235 /* default parameters (pointer to structure) */
236 /* */
237 /* error: */
238 /* return NULL */
239 /*----------------------------------------------------------------------*/
240 {
241 struct unur_par *par;
242
243 /* check arguments */
244 _unur_check_NULL( GENTYPE,distr,NULL );
245
246 /* check distribution */
247 if (distr->type != UNUR_DISTR_CVEC) {
248 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,"");
249 return NULL; }
250 COOKIE_CHECK(distr,CK_DISTR_CVEC,NULL);
251
252 if (!(distr->set & UNUR_DISTR_SET_RANKCORR)) {
253 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"rank correlation matrix");
254 return NULL; }
255
256 if (!(distr->set & UNUR_DISTR_SET_MARGINAL)) {
257 _unur_error(GENTYPE,UNUR_ERR_DISTR_REQUIRED,"marginals");
258 return NULL; }
259
260 /* allocate structure */
261 par = _unur_par_new( sizeof(struct unur_norta_par) );
262 COOKIE_SET(par,CK_NORTA_PAR);
263
264 /* copy input */
265 par->distr = distr; /* pointer to distribution object */
266
267 /* set default values */
268 par->method = UNUR_METH_NORTA ; /* method */
269 par->variant = 0u; /* default variant */
270 par->set = 0u; /* inidicate default parameters */
271 par->urng = unur_get_default_urng(); /* use default urng */
272 par->urng_aux = NULL; /* no auxilliary URNG required */
273
274 par->debug = _unur_default_debugflag; /* set default debugging flags */
275
276 /* routine for starting generator */
277 par->init = _unur_norta_init;
278
279 return par;
280
281 } /* end of unur_norta_new() */
282
283
284 /*****************************************************************************/
285 /** Private **/
286 /*****************************************************************************/
287
288 struct unur_gen *
_unur_norta_init(struct unur_par * par)289 _unur_norta_init( struct unur_par *par )
290 /*----------------------------------------------------------------------*/
291 /* initialize new generator */
292 /* */
293 /* parameters: */
294 /* par ... pointer to paramters for building generator object */
295 /* */
296 /* return: */
297 /* pointer to generator object */
298 /* */
299 /* error: */
300 /* return NULL */
301 /*----------------------------------------------------------------------*/
302 {
303 struct unur_gen *gen;
304
305 /* check arguments */
306 _unur_check_NULL( GENTYPE,par,NULL );
307
308 /* check input */
309 if ( par->method != UNUR_METH_NORTA ) {
310 _unur_error(GENTYPE,UNUR_ERR_PAR_INVALID,"");
311 return NULL; }
312 COOKIE_CHECK(par,CK_NORTA_PAR,NULL);
313
314 /* create a new empty generator object */
315 gen = _unur_norta_create(par);
316 _unur_par_free(par);
317 if (!gen) return NULL;
318
319 /* check for truncated domain */
320 if ( gen->distr->set & UNUR_DISTR_SET_DOMAINBOUNDED ) {
321 if ( DISTR.domainrect == NULL ) {
322 /* cannot handle non-rectangular domain */
323 _unur_error(gen->genid,UNUR_ERR_GEN_CONDITION,"cannot handle non-rectangular domain");
324 _unur_norta_free(gen); return NULL;
325 }
326 else { /* rectangular domain */
327 if (_unur_distr_cvec_marginals_are_equal(DISTR.marginals, GEN->dim)) {
328 /* we have to handle these equal marginal distributions independently */
329 if ( _unur_distr_cvec_duplicate_firstmarginal(gen->distr) != UNUR_SUCCESS ) {
330 _unur_norta_free(gen); return NULL;
331 }
332 }
333 }
334 }
335
336 #ifdef UNUR_ENABLE_LOGGING
337 /* write info into LOG file */
338 if (gen->debug) _unur_norta_debug_init(gen);
339 #endif
340
341 /* compute parameters for NORTU (normal to uniform) */
342 if (_unur_norta_nortu_setup(gen) != UNUR_SUCCESS) {
343 _unur_norta_free(gen); return NULL;
344 }
345
346 /* distribution object for standard normal distribution */
347 GEN->normaldistr = unur_distr_normal(NULL,0);
348
349 if (gen->distr->id != UNUR_DISTR_COPULA) {
350 /* we have to initialize generator for marginal distributions */
351
352 if (_unur_distr_cvec_marginals_are_equal(DISTR.marginals, GEN->dim)) {
353 /* we can use the same generator object for all marginal distribuitons */
354 struct unur_gen *marginalgen = _unur_norta_make_marginalgen( gen, DISTR.marginals[0] );
355 if (marginalgen)
356 GEN->marginalgen_list = _unur_gen_list_set(marginalgen,GEN->dim);
357 }
358
359 else {
360 /* we need different generator objects for all marginal distribuitons */
361 int i,j;
362 int failed = FALSE;
363 struct unur_gen **marginalgens = _unur_xmalloc( GEN->dim * sizeof(struct unur_gen*) );
364
365 /* check for truncated domain */
366 if ( gen->distr->set & UNUR_DISTR_SET_DOMAINBOUNDED ) {
367 /* set domains for each marginal distribution */
368 for (i=0; i<GEN->dim && !failed; i++) {
369 if ( (unur_distr_cont_set_domain(DISTR.marginals[i],
370 DISTR.domainrect[2*i], DISTR.domainrect[2*i+1]))
371 != UNUR_SUCCESS) {
372 failed = TRUE; break;
373 }
374 }
375 }
376
377 /* create generator for each marginal distribution */
378 for (i=0; i<GEN->dim && !failed; i++) {
379 marginalgens[i] = _unur_norta_make_marginalgen( gen, DISTR.marginals[i] );
380 if (marginalgens[i]==NULL) {
381 failed=TRUE; break;
382 }
383 }
384
385 /* check creating of generators */
386 if (failed) {
387 for (j=0; j<i; j++) _unur_free(marginalgens[j]);
388 free (marginalgens);
389 }
390 else
391 GEN->marginalgen_list = marginalgens;
392 }
393
394 /* verify initialization of marginal generators */
395 if (GEN->marginalgen_list == NULL) {
396 _unur_error(gen->genid,UNUR_ERR_GENERIC,"init of marginal generators failed");
397 _unur_norta_free(gen);
398 return NULL;
399 }
400 }
401
402 /* o.k. */
403 return gen;
404
405 } /* end of _unur_norta_init() */
406
407 /*---------------------------------------------------------------------------*/
408
409 static struct unur_gen *
_unur_norta_create(struct unur_par * par)410 _unur_norta_create( struct unur_par *par )
411 /*----------------------------------------------------------------------*/
412 /* allocate memory for generator */
413 /* */
414 /* parameters: */
415 /* par ... pointer to parameter for building generator object */
416 /* */
417 /* return: */
418 /* pointer to (empty) generator object with default settings */
419 /* */
420 /* error: */
421 /* return NULL */
422 /*----------------------------------------------------------------------*/
423 {
424 struct unur_gen *gen;
425
426 /* check arguments */
427 CHECK_NULL(par,NULL); COOKIE_CHECK(par,CK_NORTA_PAR,NULL);
428
429 /* create new generic generator object */
430 gen = _unur_generic_create( par, sizeof(struct unur_norta_gen) );
431
432 /* magic cookies */
433 COOKIE_SET(gen,CK_NORTA_GEN);
434
435 /* set generator identifier */
436 gen->genid = _unur_set_genid(GENTYPE);
437
438 /* routines for sampling and destroying generator */
439 SAMPLE = _unur_norta_getSAMPLE(gen);
440 gen->destroy = _unur_norta_free;
441 gen->clone = _unur_norta_clone;
442
443 /* dimension of distribution */
444 GEN->dim = gen->distr->dim;
445
446 /* allocate array for auxiliary copula */
447 GEN->copula = _unur_xmalloc(sizeof(double)*GEN->dim);
448
449 /* initialize pointer */
450 MNORMAL = NULL;
451 GEN->normaldistr = NULL;
452 GEN->marginalgen_list = NULL;
453
454 #ifdef UNUR_ENABLE_INFO
455 /* set function for creating info string */
456 gen->info = _unur_norta_info;
457 #endif
458
459 /* return pointer to (almost empty) generator object */
460 return gen;
461
462 } /* end of _unur_norta_create() */
463
464 /*---------------------------------------------------------------------------*/
465
466 struct unur_gen *
_unur_norta_clone(const struct unur_gen * gen)467 _unur_norta_clone( const struct unur_gen *gen )
468 /*----------------------------------------------------------------------*/
469 /* copy (clone) generator object */
470 /* */
471 /* parameters: */
472 /* gen ... pointer to generator object */
473 /* */
474 /* return: */
475 /* pointer to clone of generator object */
476 /* */
477 /* error: */
478 /* return NULL */
479 /*----------------------------------------------------------------------*/
480 {
481 #define CLONE ((struct unur_norta_gen*)clone->datap)
482
483 struct unur_gen *clone;
484
485 /* check arguments */
486 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_NORTA_GEN,NULL);
487
488 /* create generic clone */
489 clone = _unur_generic_clone( gen, GENTYPE );
490
491 /* allocate array for auxiliary copula */
492 CLONE->copula = _unur_xmalloc(sizeof(double)*GEN->dim);
493
494 /* clone marginal distribution */
495 CLONE->normaldistr = _unur_distr_clone(GEN->normaldistr);
496
497 /* marginal generators */
498 if (GEN->marginalgen_list)
499 CLONE->marginalgen_list = _unur_gen_list_clone( GEN->marginalgen_list, GEN->dim );
500
501 return clone;
502
503 #undef CLONE
504 } /* end of _unur_norta_clone() */
505
506 /*---------------------------------------------------------------------------*/
507
508 void
_unur_norta_free(struct unur_gen * gen)509 _unur_norta_free( struct unur_gen *gen )
510 /*----------------------------------------------------------------------*/
511 /* deallocate generator object */
512 /* */
513 /* parameters: */
514 /* gen ... pointer to generator object */
515 /*----------------------------------------------------------------------*/
516 {
517 /* check arguments */
518 if( !gen ) /* nothing to do */
519 return;
520
521 /* check input */
522 if ( gen->method != UNUR_METH_NORTA ) {
523 _unur_warning(gen->genid,UNUR_ERR_GEN_INVALID,"");
524 return; }
525 COOKIE_CHECK(gen,CK_NORTA_GEN,RETURN_VOID);
526
527 /* free auxiliary arrays */
528 if (GEN->copula) free (GEN->copula);
529
530 /* free normal distribution object */
531 if (GEN->normaldistr) _unur_distr_free (GEN->normaldistr);
532
533 /* free marginal generators */
534 if (GEN->marginalgen_list)
535 _unur_gen_list_free( GEN->marginalgen_list, GEN->dim);
536
537 /* we cannot use this generator object any more */
538 SAMPLE = NULL; /* make sure to show up a programming error */
539
540 _unur_generic_free(gen);
541
542 } /* end of _unur_norta_free() */
543
544 /*****************************************************************************/
545
546 int
_unur_norta_sample_cvec(struct unur_gen * gen,double * vec)547 _unur_norta_sample_cvec( struct unur_gen *gen, double *vec )
548 /*----------------------------------------------------------------------*/
549 /* sample from generator */
550 /* */
551 /* parameters: */
552 /* gen ... pointer to generator object */
553 /* vec ... random vector (result) */
554 /*----------------------------------------------------------------------*/
555 {
556 #define idx(a,b) ((a)*GEN->dim+(b))
557 int j;
558 double *u;
559
560 /* check arguments */
561 CHECK_NULL(gen,UNUR_ERR_NULL);
562 COOKIE_CHECK(gen,CK_NORTA_GEN,UNUR_ERR_COOKIE);
563
564 /* pointer to auxiliary array of uniforms */
565 u = GEN->copula;
566
567 /* sample from multinormal distribution */
568 _unur_sample_vec(MNORMAL,u);
569
570 /* make copula */
571 for (j=0; j<GEN->dim; j++)
572 vec[j] = unur_distr_cont_eval_cdf( u[j], GEN->normaldistr );
573
574 if (gen->distr->id == UNUR_DISTR_COPULA)
575 /* we want to have a normal copula --> just return data */
576 return UNUR_SUCCESS;
577
578 /* else non-uniform marginals */
579 for (j=0; j<GEN->dim; j++) {
580 /* call marginal generator (using this U-value) */
581 vec[j] = unur_quantile(GEN->marginalgen_list[j], vec[j]);
582 }
583
584 return UNUR_SUCCESS;
585
586 #undef idx
587 } /* end of _unur_norta_sample_cvec() */
588
589
590 /*****************************************************************************/
591 /** Auxilliary Routines **/
592 /*****************************************************************************/
593
594 int
_unur_norta_nortu_setup(struct unur_gen * gen)595 _unur_norta_nortu_setup( struct unur_gen *gen )
596 /*----------------------------------------------------------------------*/
597 /* Compute parameter for NORTU (Normal to uniform) */
598 /* */
599 /* parameters: */
600 /* gen ... pointer to generator object */
601 /* */
602 /* return: */
603 /* UNUR_SUCCESS if computed successfully */
604 /* error code otherwise */
605 /*----------------------------------------------------------------------*/
606 {
607 #define idx(a,b) ((a)*dim+(b))
608
609 int dim = GEN->dim; /* dimension of distribution */
610 double *sigma_y; /* correlation matrix for corresponding normal distr. */
611 double *eigenvalues; /* eigenvalues of sigma_y */
612 double *eigenvectors; /* eigenvectors of sigma_y */
613 int eigenvalues_positive; /* boolean indicating whether all eigenvalues are
614 strictly positive */
615 struct unur_distr *mn_distr; /* multinormal distribution */
616 struct unur_gen *mn_gen; /* generator for multinormal distribution */
617 int i,j;
618
619 /* setup correlation matrix for corresponding normal distribution */
620 sigma_y = _unur_xmalloc(dim * dim * sizeof(double));
621 for(i=0; i<dim; i++) {
622 /* left lower part: make matrix symmetric */
623 for(j=0; j<i; j++)
624 sigma_y[idx(i,j)] = sigma_y[idx(j,i)];
625 /* diagonal */
626 sigma_y[idx(i,i)] = 1.;
627 /* right upper part */
628 for(j=i+1; j<dim; j++)
629 sigma_y[idx(i,j)] = 2.*sin(DISTR.rankcorr[idx(i,j)]*(M_PI/6.));
630 }
631
632 #ifdef UNUR_ENABLE_LOGGING
633 /* write info into LOG file */
634 if (gen->debug & NORTA_DEBUG_SIGMA_Y)
635 _unur_norta_debug_sigma_y( gen, sigma_y, "NORTU setup:" );
636 #endif
637
638 /* compute eigenvectors and eigenvalues of sigma_y */
639 eigenvalues = _unur_xmalloc(dim * sizeof(double));
640 eigenvectors = _unur_xmalloc(dim * dim * sizeof(double));
641 if (_unur_matrix_eigensystem(dim, sigma_y, eigenvalues, eigenvectors) != UNUR_SUCCESS) {
642 _unur_error(GENTYPE,UNUR_ERR_GEN_DATA,"cannot compute eigenvalues for given sigma_y");
643 free(sigma_y); free(eigenvalues); free(eigenvectors);
644 return UNUR_ERR_GEN_DATA;
645 }
646
647 #ifdef UNUR_ENABLE_LOGGING
648 if (gen->debug & NORTA_DEBUG_SIGMA_Y)
649 _unur_norta_debug_eigensystem( gen, eigenvalues, eigenvectors );
650 #endif
651
652 /* check if all eigenvalues are positive */
653 /* otherwise set to small values close to 0 */
654 eigenvalues_positive = TRUE;
655 for(i=0; i<dim; i++)
656 if(eigenvalues[i] < UNUR_NORTA_MIN_EIGENVALUE) {
657 eigenvalues[i] = UNUR_NORTA_MIN_EIGENVALUE;
658 eigenvalues_positive = FALSE;
659 }
660
661 /* make corrected correlation matrix */
662
663 if (!eigenvalues_positive) {
664 _unur_matrix_transform_diagonal(dim,eigenvectors,eigenvalues,sigma_y);
665 _unur_norta_make_correlationmatrix(dim,sigma_y);
666 _unur_warning(GENTYPE,UNUR_ERR_GEN_DATA,
667 "sigma_y not positive definite -> corrected matrix");
668 #ifdef UNUR_ENABLE_LOGGING
669 if (gen->debug & NORTA_DEBUG_SIGMA_Y)
670 _unur_norta_debug_sigma_y( gen, sigma_y, "\tEigenvalues < 0 --> correction required" );
671 #endif
672 }
673
674 /* clear working arrays */
675 free(eigenvalues);
676 free(eigenvectors);
677
678 /* make generator for multinormal distribution */
679 mn_distr = unur_distr_multinormal(dim, NULL, sigma_y);
680 mn_gen = NULL;
681 if (mn_distr) {
682 mn_gen = unur_init(unur_mvstd_new(mn_distr));
683 _unur_distr_free(mn_distr);
684 }
685 if (mn_gen == NULL) {
686 _unur_error(GENTYPE,UNUR_ERR_GEN_DATA,"(corrected) sigma_y not positive definit");
687 free(sigma_y);
688 return UNUR_ERR_GEN_DATA;
689 }
690 MNORMAL = mn_gen;
691 /* need same uniform random number generator as NORTA generator */
692 MNORMAL->urng = gen->urng;
693 /* copy debugging flags */
694 MNORMAL->debug = gen->debug;
695
696 #ifdef UNUR_ENABLE_LOGGING
697 if (gen->debug & NORTA_DEBUG_SIGMA_Y)
698 _unur_norta_debug_nmgenerator( gen );
699 #endif
700
701 /* clear working arrays */
702 free(sigma_y);
703
704 return UNUR_SUCCESS;
705
706 #undef idx
707 } /* end of _unur_norta_nortu_setup() */
708
709 /*---------------------------------------------------------------------------*/
710
711 int
_unur_norta_make_correlationmatrix(int dim,double * M)712 _unur_norta_make_correlationmatrix( int dim, double *M)
713 /*----------------------------------------------------------------------*/
714 /* make correlation matrix by transforming symmetric positive definit */
715 /* matrix. */
716 /* */
717 /* There is no checking whether M fulfills the conditions! */
718 /* */
719 /* parameters: */
720 /* dim ... dimension of matrix M */
721 /* M ... symmetric square matrix */
722 /*----------------------------------------------------------------------*/
723 {
724 #define idx(a,b) ((a)*dim+(b))
725
726 int i,j;
727
728 /* diagonal is used to store the roots of the diagonal elements */
729 for (i=0; i<dim; i++)
730 M[idx(i,i)] = sqrt(M[idx(i,i)]);
731
732 for (i=0; i<dim; i++)
733 for (j=i; j<dim; j++)
734 if(i!=j) {
735 M[idx(i,j)] /= M[idx(i,i)] * M[idx(j,j)];
736 M[idx(j,i)] = M[idx(i,j)];
737 }
738
739 /* the diagonal elements are set to 1. */
740 for (i=0; i<dim; i++)
741 M[idx(i,i)] = 1.;
742
743 return UNUR_SUCCESS;
744 #undef idx
745 } /* end of _unur_norta_make_correlationmatrix() */
746
747 /*---------------------------------------------------------------------------*/
748
749 struct unur_gen *
_unur_norta_make_marginalgen(const struct unur_gen * gen,const struct unur_distr * marginal)750 _unur_norta_make_marginalgen( const struct unur_gen *gen,
751 const struct unur_distr *marginal )
752 /*----------------------------------------------------------------------*/
753 /* make generator for marginal distribution */
754 /* */
755 /* parameters: */
756 /* gen ... pointer to multivariate generator object */
757 /* marginal ... distribution object for marginal distribution */
758 /* */
759 /* return: */
760 /* generator for marginal distribution */
761 /* */
762 /* error: */
763 /* return NULL */
764 /*----------------------------------------------------------------------*/
765 {
766 struct unur_gen *marginalgen;
767 struct unur_par *par;
768
769 /* check arguments */
770 CHECK_NULL(gen,NULL); COOKIE_CHECK(gen,CK_NORTA_GEN,NULL);
771 CHECK_NULL(marginal,NULL);
772
773 /* check distribution */
774 if (marginal->type != UNUR_DISTR_CONT) {
775 _unur_error(GENTYPE,UNUR_ERR_DISTR_INVALID,""); return NULL; }
776 COOKIE_CHECK(marginal,CK_DISTR_CONT,NULL);
777
778 /* try PINV, CSTD+Inversion, HINV, and NINV */
779 do {
780 /* PINV (Polynomial interpolation based INVersion of CDF) */
781 par = unur_pinv_new( marginal );
782 if ( (marginalgen = _unur_init(par)) != NULL )
783 break;
784
785 /* CSTD + inversion */
786 par = unur_cstd_new( marginal );
787 if (unur_cstd_set_variant( par, UNUR_STDGEN_INVERSION)==UNUR_SUCCESS) {
788 marginalgen = _unur_init(par);
789 break;
790 }
791 else {
792 _unur_par_free(par);
793 }
794
795 /* HINV (Hermite interpolation based INVersion of CDF) */
796 par = unur_hinv_new( marginal );
797 if ( (marginalgen = _unur_init(par)) != NULL )
798 break;
799
800 /* NINV (Numerical inversion with regula falsi */
801 par = unur_ninv_new( marginal );
802 unur_ninv_set_table( par, 100 );
803 if ( (marginalgen = _unur_init(par)) != NULL )
804 break;
805
806 /* no inversion method avaiblable */
807 _unur_error(gen->genid,UNUR_ERR_DISTR_REQUIRED,"data for (numerical) inversion of marginal missing");
808 return NULL;
809
810 } while (1);
811
812 /* set debugging flags */
813 marginalgen->debug = gen->debug;
814
815 /* return generator object */
816 return marginalgen;
817
818 } /* end of _unur_norta_make_marginalgen() */
819
820
821 /*****************************************************************************/
822 /** Debugging utilities **/
823 /*****************************************************************************/
824
825 /*---------------------------------------------------------------------------*/
826 #ifdef UNUR_ENABLE_LOGGING
827 /*---------------------------------------------------------------------------*/
828
829 void
_unur_norta_debug_init(const struct unur_gen * gen)830 _unur_norta_debug_init( const struct unur_gen *gen )
831 /*----------------------------------------------------------------------*/
832 /* write info about generator into LOG file */
833 /* */
834 /* parameters: */
835 /* gen ... pointer to generator object */
836 /*----------------------------------------------------------------------*/
837 {
838 /* int i; */
839 FILE *LOG;
840
841 /* check arguments */
842 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_NORTA_GEN,RETURN_VOID);
843
844 LOG = unur_get_stream();
845
846 fprintf(LOG,"%s:\n",gen->genid);
847 fprintf(LOG,"%s: type = continuous multivariate random variates\n",gen->genid);
848 fprintf(LOG,"%s: method = NORTA (Vector Matrix Transformation)\n",gen->genid);
849 fprintf(LOG,"%s:\n",gen->genid);
850
851 _unur_distr_cvec_debug( gen->distr, gen->genid );
852
853 } /* end of _unur_norta_debug_init() */
854
855 /*---------------------------------------------------------------------------*/
856
857 void
_unur_norta_debug_sigma_y(const struct unur_gen * gen,const double * sigma_y,const char * comment)858 _unur_norta_debug_sigma_y( const struct unur_gen *gen,
859 const double *sigma_y,
860 const char *comment )
861 /*----------------------------------------------------------------------*/
862 /* print sigma_y of corresponding normal distribution. */
863 /* */
864 /* parameters: */
865 /* gen ... pointer to generator object */
866 /* sigma_y ... pointer to correlation matrix */
867 /* comment ... additional string printed */
868 /*----------------------------------------------------------------------*/
869 {
870 FILE *LOG;
871
872 /* check arguments */
873 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_NORTA_GEN,RETURN_VOID);
874
875 LOG = unur_get_stream();
876
877 fprintf(LOG,"%s: %s\n",gen->genid,comment);
878 fprintf(LOG,"%s:\n",gen->genid);
879 _unur_matrix_print_matrix( GEN->dim, sigma_y, "\tsigma_y =",
880 LOG, gen->genid, "\t ");
881
882 } /* end of _unur_norta_debug_sigma_y() */
883
884 /*---------------------------------------------------------------------------*/
885
886 void
_unur_norta_debug_eigensystem(const struct unur_gen * gen,const double * eigenvalues,const double * eigenvectors)887 _unur_norta_debug_eigensystem( const struct unur_gen *gen,
888 const double *eigenvalues,
889 const double *eigenvectors )
890 /*----------------------------------------------------------------------*/
891 /* print eigensystem of sigma_y. */
892 /* */
893 /* parameters: */
894 /* gen ... pointer to generator object */
895 /* eigenvalues ... eigenvalues */
896 /* eigenvectors ... eigenvalues */
897 /*----------------------------------------------------------------------*/
898 {
899 FILE *LOG;
900
901 /* check arguments */
902 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_NORTA_GEN,RETURN_VOID);
903
904 LOG = unur_get_stream();
905
906 _unur_matrix_print_vector( GEN->dim, eigenvalues,
907 "\teigenvalues of sigma_y =",
908 LOG, gen->genid, "\t ");
909 _unur_matrix_print_matrix( GEN->dim, eigenvectors,
910 "\teigenvectors of sigma_y [rows] =",
911 LOG, gen->genid, "\t ");
912
913 } /* end of _unur_norta_debug_eigensystem() */
914
915 /*---------------------------------------------------------------------------*/
916
917 void
_unur_norta_debug_nmgenerator(const struct unur_gen * gen)918 _unur_norta_debug_nmgenerator( const struct unur_gen *gen )
919 /*----------------------------------------------------------------------*/
920 /* print genid of multinormal generator. */
921 /* */
922 /* parameters: */
923 /* gen ... pointer to generator object */
924 /*----------------------------------------------------------------------*/
925 {
926 FILE *LOG;
927
928 /* check arguments */
929 CHECK_NULL(gen,RETURN_VOID); COOKIE_CHECK(gen,CK_NORTA_GEN,RETURN_VOID);
930
931 LOG = unur_get_stream();
932
933 fprintf(LOG,"%s: generator for multinormal auxiliary distribution = %s\n", gen->genid,
934 MNORMAL->genid );
935 fprintf(LOG,"%s:\n",gen->genid);
936
937 } /* end of _unur_norta_debug_nmgenerator() */
938
939 /*---------------------------------------------------------------------------*/
940 #endif /* UNUR_ENABLE_LOGGING */
941 /*---------------------------------------------------------------------------*/
942
943
944 /*---------------------------------------------------------------------------*/
945 #ifdef UNUR_ENABLE_INFO
946 /*---------------------------------------------------------------------------*/
947
948 void
_unur_norta_info(struct unur_gen * gen,int help)949 _unur_norta_info( struct unur_gen *gen, int help )
950 /*----------------------------------------------------------------------*/
951 /* create character string that contains information about the */
952 /* given generator object. */
953 /* */
954 /* parameters: */
955 /* gen ... pointer to generator object */
956 /* help ... whether to print additional comments */
957 /*----------------------------------------------------------------------*/
958 {
959 struct unur_string *info = gen->infostr;
960 struct unur_distr *distr = gen->distr;
961 int i;
962
963 /* generator ID */
964 _unur_string_append(info,"generator ID: %s\n\n", gen->genid);
965
966 /* distribution */
967 _unur_string_append(info,"distribution:\n");
968 _unur_distr_info_typename(gen);
969 _unur_string_append(info," dimension = %d\n",GEN->dim);
970 _unur_string_append(info," functions = MARGINAL distributions\n");
971
972 _unur_string_append(info," marginals =");
973 for (i=0; i<distr->dim; i++)
974 _unur_string_append(info," %s", distr->data.cvec.marginals[i]->name);
975 _unur_string_append(info,"\n\n");
976
977 /* if (help) { */
978 /* _unur_string_append(info,"\n"); */
979 /* } */
980
981 /* method */
982 _unur_string_append(info,"method: NORTA (NORmal To Anything)\n");
983 _unur_string_append(info,"\n");
984
985 /* performance */
986 /* _unur_string_append(info,"performance characteristics:\n"); */
987 /* _unur_string_append(info,"\n"); */
988
989 /* parameters */
990 if (help) {
991 _unur_string_append(info,"parameters: none\n");
992 _unur_string_append(info,"\n");
993 }
994
995 /* Hints */
996 /* if (help) { */
997 /* _unur_string_append(info,"\n"); */
998 /* } */
999
1000 } /* end of _unur_norta_info() */
1001
1002 /*---------------------------------------------------------------------------*/
1003 #endif /* end UNUR_ENABLE_INFO */
1004 /*---------------------------------------------------------------------------*/
1005