1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: cvec.c *
8 * *
9 * manipulate multivariate continuous distribution objects *
10 * *
11 * return: *
12 * UNUR_SUCCESS ... on success *
13 * error code ... on error *
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 /*---------------------------------------------------------------------------*/
38
39 #include <unur_source.h>
40 #include <distributions/unur_stddistr.h>
41 #include "distr_source.h"
42 #include "distr.h"
43 #include "cvec.h"
44 #include <utils/matrix_source.h>
45 #include <stdarg.h>
46
47 /*---------------------------------------------------------------------------*/
48 /* copy and free marginal distributions */
49
50 static struct unur_distr **_unur_distr_cvec_marginals_clone ( struct unur_distr **marginals, int dim );
51 static void _unur_distr_cvec_marginals_free ( struct unur_distr **marginals, int dim );
52
53 static void _unur_distr_cvec_free( struct unur_distr *distr );
54
55 /*---------------------------------------------------------------------------*/
56
57 #define DISTR distr->data.cvec
58
59 /*---------------------------------------------------------------------------*/
60
61 /*****************************************************************************/
62 /** **/
63 /** mulitvariate continuous distributions **/
64 /** **/
65 /*****************************************************************************/
66
67 /*---------------------------------------------------------------------------*/
68
69 struct unur_distr *
unur_distr_cvec_new(int dim)70 unur_distr_cvec_new( int dim )
71 /*----------------------------------------------------------------------*/
72 /* create a new (empty) distribution object */
73 /* type: multivariate continuous with given PDF */
74 /* */
75 /* parameters: */
76 /* dim ... number of components of random vector (dimension) */
77 /* */
78 /* return: */
79 /* pointer to distribution object */
80 /* */
81 /* error: */
82 /* return NULL */
83 /*----------------------------------------------------------------------*/
84 {
85 register struct unur_distr *distr;
86 int i;
87
88 /* check dimension for new parameter for distribution */
89 if (dim < 1) {
90 _unur_error(NULL,UNUR_ERR_DISTR_SET,"dimension < 1");
91 return NULL;
92 }
93
94 /* get empty distribution object */
95 distr = _unur_distr_generic_new();
96 if (!distr) return NULL;
97
98 /* set magic cookie */
99 COOKIE_SET(distr,CK_DISTR_CVEC);
100
101 /* set type of distribution */
102 distr->type = UNUR_DISTR_CVEC;
103
104 /* set id to generic distribution */
105 distr->id = UNUR_DISTR_GENERIC;
106
107 /* dimension of random vector */
108 distr->dim = dim; /* multivariant */
109
110 /* this is not a derived distribution */
111 distr->base = NULL;
112
113 /* destructor */
114 distr->destroy = _unur_distr_cvec_free;
115
116 /* clone */
117 distr->clone = _unur_distr_cvec_clone;
118
119 /* set defaults */
120 DISTR.pdf = NULL; /* pointer to PDF */
121 DISTR.dpdf = NULL; /* pointer to gradient of PDF */
122 DISTR.pdpdf = NULL; /* pointer to partial derivative of PDF */
123 DISTR.logpdf = NULL; /* pointer to logPDF */
124 DISTR.dlogpdf = NULL; /* pointer to gradient of logPDF */
125 DISTR.pdlogpdf = NULL; /* pointer to partial derivative of logPDF */
126 DISTR.domainrect = NULL; /* (rectangular) domain of distribution [default: unbounded */
127 DISTR.init = NULL; /* pointer to special init routine [default: none] */
128 DISTR.mean = NULL; /* mean vector [default: not known] */
129 DISTR.covar = NULL; /* covariance matrix [default: not known] */
130 DISTR.cholesky = NULL; /* cholesky factor of cov. matrix [default: not computed] */
131 DISTR.covar_inv = NULL; /* inverse covariance matrix [default: not computed] */
132 DISTR.rankcorr = NULL; /* rank correlation [default: not known] */
133 DISTR.rk_cholesky = NULL; /* cholesky factor of rank correlation [default: not computed] */
134 DISTR.marginals = NULL; /* array of pointers to marginal distributions */
135 DISTR.upd_mode = NULL; /* funct for computing mode */
136 DISTR.upd_volume = NULL; /* funct for computing volume */
137
138 #ifdef USE_DEPRECATED_CODE
139 DISTR.stdmarginals = NULL; /* array of pointers to standardized marginal distributions */
140 #endif
141
142 /* initialize parameters of the p.d.f. */
143 DISTR.n_params = 0; /* number of parameters of the pdf */
144 for (i=0; i<UNUR_DISTR_MAXPARAMS; i++)
145 DISTR.params[i] = 0.;
146
147 /* initialize parameter vectors of the PDF */
148 for (i=0; i<UNUR_DISTR_MAXPARAMS; i++) {
149 DISTR.n_param_vec[i] = 0;
150 DISTR.param_vecs[i] = NULL;
151 }
152
153 DISTR.norm_constant = 1.; /* (log of) normalization constant for PDF
154 (initialized to avoid accidently floating
155 point exception */
156
157 DISTR.mode = NULL; /* location of mode (default: not known) */
158 DISTR.center = NULL; /* location of center (default: not given)*/
159 DISTR.volume = INFINITY; /* area below PDF (default: not known) */
160
161 /* return pointer to object */
162 return distr;
163
164 } /* end of unur_distr_cvec_new() */
165
166 /*---------------------------------------------------------------------------*/
167
168 struct unur_distr *
_unur_distr_cvec_clone(const struct unur_distr * distr)169 _unur_distr_cvec_clone( const struct unur_distr *distr )
170 /*----------------------------------------------------------------------*/
171 /* copy (clone) distribution object */
172 /* */
173 /* parameters: */
174 /* distr ... pointer to source distribution object */
175 /* */
176 /* return: */
177 /* pointer to clone of distribution object */
178 /*----------------------------------------------------------------------*/
179 {
180 #define CLONE clone->data.cvec
181
182 struct unur_distr *clone;
183 int i;
184
185 /* check arguments */
186 _unur_check_NULL( NULL, distr, NULL );
187 _unur_check_distr_object( distr, CVEC, NULL );
188
189 /* allocate memory */
190 clone = _unur_xmalloc( sizeof(struct unur_distr) );
191
192 /* copy distribution object into clone */
193 memcpy( clone, distr, sizeof( struct unur_distr ) );
194
195 /* copy data about distribution */
196 if (DISTR.domainrect) {
197 CLONE.domainrect = _unur_xmalloc( 2 * distr->dim * sizeof(double) );
198 memcpy( CLONE.domainrect, DISTR.domainrect, 2 * distr->dim * sizeof(double) );
199 }
200
201 if (DISTR.mean) {
202 CLONE.mean = _unur_xmalloc( distr->dim * sizeof(double) );
203 memcpy( CLONE.mean, DISTR.mean, distr->dim * sizeof(double) );
204 }
205
206 if (DISTR.covar) {
207 CLONE.covar = _unur_xmalloc( distr->dim * distr->dim * sizeof(double) );
208 memcpy( CLONE.covar, DISTR.covar, distr->dim * distr->dim * sizeof(double) );
209 }
210
211 if (DISTR.cholesky) {
212 CLONE.cholesky = _unur_xmalloc( distr->dim * distr->dim * sizeof(double) );
213 memcpy( CLONE.cholesky, DISTR.cholesky, distr->dim * distr->dim * sizeof(double) );
214 }
215
216 if (DISTR.covar_inv) {
217 CLONE.covar_inv = _unur_xmalloc( distr->dim * distr->dim * sizeof(double) );
218 memcpy( CLONE.covar_inv, DISTR.covar_inv, distr->dim * distr->dim * sizeof(double) );
219 }
220
221 if (DISTR.rankcorr) {
222 CLONE.rankcorr = _unur_xmalloc( distr->dim * distr->dim * sizeof(double) );
223 memcpy( CLONE.rankcorr, DISTR.rankcorr, distr->dim * distr->dim * sizeof(double) );
224 }
225
226 if (DISTR.rk_cholesky) {
227 CLONE.rk_cholesky = _unur_xmalloc( distr->dim * distr->dim * sizeof(double) );
228 memcpy( CLONE.rk_cholesky, DISTR.rk_cholesky, distr->dim * distr->dim * sizeof(double) );
229 }
230
231 if (DISTR.mode) {
232 CLONE.mode = _unur_xmalloc( distr->dim * sizeof(double) );
233 memcpy( CLONE.mode, DISTR.mode, distr->dim * sizeof(double) );
234 }
235
236 if (DISTR.center) {
237 CLONE.center = _unur_xmalloc( distr->dim * sizeof(double) );
238 memcpy( CLONE.center, DISTR.center, distr->dim * sizeof(double) );
239 }
240
241 if (DISTR.marginals)
242 CLONE.marginals = _unur_distr_cvec_marginals_clone( DISTR.marginals, distr->dim );
243
244 #ifdef USE_DEPRECATED_CODE
245 if (DISTR.stdmarginals)
246 CLONE.stdmarginals = _unur_distr_cvec_marginals_clone( DISTR.stdmarginals, distr->dim );
247 #endif
248
249 /* clone of scalar parameters */
250 CLONE.n_params = DISTR.n_params;
251 for (i=0; i<UNUR_DISTR_MAXPARAMS; i++) {
252 CLONE.params[i] = DISTR.params[i];
253 }
254
255 /* clone of parameter arrays */
256 for (i=0; i<UNUR_DISTR_MAXPARAMS; i++) {
257 CLONE.n_param_vec[i] = DISTR.n_param_vec[i];
258 if (DISTR.param_vecs[i]) {
259 CLONE.param_vecs[i] = _unur_xmalloc( DISTR.n_param_vec[i] * sizeof(double) );
260 memcpy( CLONE.param_vecs[i], DISTR.param_vecs[i], DISTR.n_param_vec[i] * sizeof(double) );
261 }
262 }
263
264 /* copy user name for distribution */
265 if (distr->name_str) {
266 size_t len = strlen(distr->name_str) + 1;
267 clone->name_str = _unur_xmalloc(len);
268 memcpy( clone->name_str, distr->name_str, len );
269 clone->name = clone->name_str;
270 }
271
272 return clone;
273
274 #undef CLONE
275 } /* end of _unur_distr_cvec_clone() */
276
277 /*---------------------------------------------------------------------------*/
278
279 void
_unur_distr_cvec_free(struct unur_distr * distr)280 _unur_distr_cvec_free( struct unur_distr *distr )
281 /*----------------------------------------------------------------------*/
282 /* free distribution object */
283 /* */
284 /* parameters: */
285 /* distr ... pointer to distribution object */
286 /*----------------------------------------------------------------------*/
287 {
288 int i;
289
290 /* check arguments */
291 if( distr == NULL ) /* nothing to do */
292 return;
293
294 COOKIE_CHECK(distr,CK_DISTR_CVEC,RETURN_VOID);
295
296 for (i=0; i<UNUR_DISTR_MAXPARAMS; i++)
297 if (DISTR.param_vecs[i]) free( DISTR.param_vecs[i] );
298
299 if (DISTR.domainrect) free(DISTR.domainrect);
300 if (DISTR.mean) free(DISTR.mean);
301 if (DISTR.covar) free(DISTR.covar);
302 if (DISTR.covar_inv) free(DISTR.covar_inv);
303 if (DISTR.cholesky) free(DISTR.cholesky);
304 if (DISTR.rankcorr) free(DISTR.rankcorr);
305 if (DISTR.rk_cholesky) free(DISTR.rk_cholesky);
306
307 if (DISTR.mode) free(DISTR.mode);
308 if (DISTR.center) free(DISTR.center);
309
310 if (DISTR.marginals)
311 _unur_distr_cvec_marginals_free(DISTR.marginals, distr->dim);
312
313 #ifdef USE_DEPRECATED_CODE
314 if (DISTR.stdmarginals)
315 _unur_distr_cvec_marginals_free(DISTR.stdmarginals, distr->dim);
316 #endif
317
318 /* user name for distribution */
319 if (distr->name_str) free(distr->name_str);
320
321 COOKIE_CLEAR(distr);
322 free( distr );
323
324 } /* end of unur_distr_cvec_free() */
325
326 /*---------------------------------------------------------------------------*/
327
328 int
unur_distr_cvec_set_pdf(struct unur_distr * distr,UNUR_FUNCT_CVEC * pdf)329 unur_distr_cvec_set_pdf( struct unur_distr *distr, UNUR_FUNCT_CVEC *pdf )
330 /*----------------------------------------------------------------------*/
331 /* set PDF of distribution */
332 /* */
333 /* parameters: */
334 /* distr ... pointer to distribution object */
335 /* pdf ... pointer to PDF */
336 /* */
337 /* return: */
338 /* UNUR_SUCCESS ... on success */
339 /* error code ... on error */
340 /*----------------------------------------------------------------------*/
341 {
342 /* check arguments */
343 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
344 _unur_check_NULL( distr->name, pdf, UNUR_ERR_NULL);
345 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
346
347 /* we do not allow overwriting a PDF */
348 if (DISTR.pdf != NULL || DISTR.logpdf != NULL) {
349 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of PDF not allowed");
350 return UNUR_ERR_DISTR_SET;
351 }
352
353 /* changelog */
354 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
355 /* derived parameters like mode, area, etc. might be wrong now! */
356
357 DISTR.pdf = pdf;
358 return UNUR_SUCCESS;
359
360 } /* end of unur_distr_cvec_set_pdf() */
361
362 /*---------------------------------------------------------------------------*/
363
364 int
unur_distr_cvec_set_dpdf(struct unur_distr * distr,UNUR_VFUNCT_CVEC * dpdf)365 unur_distr_cvec_set_dpdf( struct unur_distr *distr, UNUR_VFUNCT_CVEC *dpdf )
366 /*----------------------------------------------------------------------*/
367 /* set gradient of PDF of distribution */
368 /* */
369 /* parameters: */
370 /* distr ... pointer to distribution object */
371 /* dpdf ... pointer to gradient of PDF */
372 /* */
373 /* return: */
374 /* UNUR_SUCCESS ... on success */
375 /* error code ... on error */
376 /*----------------------------------------------------------------------*/
377 {
378 /* check arguments */
379 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
380 _unur_check_NULL( distr->name, dpdf, UNUR_ERR_NULL );
381 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
382
383 /* we do not allow overwriting a dPDF */
384 if (DISTR.dpdf != NULL || DISTR.dlogpdf != NULL) {
385 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of dPDF not allowed");
386 return UNUR_ERR_DISTR_SET;
387 }
388
389 /* changelog */
390 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
391 /* derived parameters like mode, area, etc. might be wrong now! */
392
393 DISTR.dpdf = dpdf;
394 return UNUR_SUCCESS;
395 } /* end of unur_distr_cvec_set_dpdf() */
396
397 /*---------------------------------------------------------------------------*/
398
399 int
unur_distr_cvec_set_pdpdf(struct unur_distr * distr,UNUR_FUNCTD_CVEC * pdpdf)400 unur_distr_cvec_set_pdpdf( struct unur_distr *distr, UNUR_FUNCTD_CVEC *pdpdf )
401 /*----------------------------------------------------------------------*/
402 /* set PDF of distribution */
403 /* */
404 /* parameters: */
405 /* distr ... pointer to distribution object */
406 /* pdpdf ... pointer to partial derivative of the PDF */
407 /* */
408 /* return: */
409 /* UNUR_SUCCESS ... on success */
410 /* error code ... on error */
411 /*----------------------------------------------------------------------*/
412 {
413 /* check arguments */
414 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
415 _unur_check_NULL( distr->name, pdpdf, UNUR_ERR_NULL);
416 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
417
418 /* we do not allow overwriting a pdPDF */
419 if (DISTR.pdpdf != NULL || DISTR.pdlogpdf != NULL) {
420 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of pdPDF not allowed");
421 return UNUR_ERR_DISTR_SET;
422 }
423
424 /* changelog */
425 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
426 /* derived parameters like mode, area, etc. might be wrong now! */
427
428 DISTR.pdpdf = pdpdf;
429 return UNUR_SUCCESS;
430
431 } /* end of unur_distr_cvec_set_pdpdf() */
432
433 /*---------------------------------------------------------------------------*/
434
435 UNUR_FUNCT_CVEC *
unur_distr_cvec_get_pdf(const struct unur_distr * distr)436 unur_distr_cvec_get_pdf( const struct unur_distr *distr )
437 /*----------------------------------------------------------------------*/
438 /* get pointer to PDF of distribution */
439 /* */
440 /* parameters: */
441 /* distr ... pointer to distribution object */
442 /* */
443 /* return: */
444 /* pointer to PDF */
445 /*----------------------------------------------------------------------*/
446 {
447 /* check arguments */
448 _unur_check_NULL( NULL, distr, NULL );
449 _unur_check_distr_object( distr, CVEC, NULL );
450
451 return DISTR.pdf;
452 } /* end of unur_distr_cvec_get_pdf() */
453
454 /*---------------------------------------------------------------------------*/
455
456 UNUR_VFUNCT_CVEC *
unur_distr_cvec_get_dpdf(const struct unur_distr * distr)457 unur_distr_cvec_get_dpdf( const struct unur_distr *distr )
458 /*----------------------------------------------------------------------*/
459 /* get pointer to gradient of PDF of distribution */
460 /* */
461 /* parameters: */
462 /* distr ... pointer to distribution object */
463 /* */
464 /* return: */
465 /* pointer to gradient of PDF */
466 /*----------------------------------------------------------------------*/
467 {
468 /* check arguments */
469 _unur_check_NULL( NULL, distr, NULL );
470 _unur_check_distr_object( distr, CVEC, NULL );
471
472 return DISTR.dpdf;
473 } /* end of unur_distr_cvec_get_dpdf() */
474
475 /*---------------------------------------------------------------------------*/
476
477 UNUR_FUNCTD_CVEC *
unur_distr_cvec_get_pdpdf(const struct unur_distr * distr)478 unur_distr_cvec_get_pdpdf( const struct unur_distr *distr )
479 /*----------------------------------------------------------------------*/
480 /* get pointer to partial derivative of PDF of distribution */
481 /* */
482 /* parameters: */
483 /* distr ... pointer to distribution object */
484 /* */
485 /* return: */
486 /* pointer to partial derivative of PDF */
487 /*----------------------------------------------------------------------*/
488 {
489 /* check arguments */
490 _unur_check_NULL( NULL, distr, NULL );
491 _unur_check_distr_object( distr, CVEC, NULL );
492
493 return DISTR.pdpdf;
494 } /* end of unur_distr_cvec_get_pdpdf() */
495
496 /*---------------------------------------------------------------------------*/
497
498 double
unur_distr_cvec_eval_pdf(const double * x,struct unur_distr * distr)499 unur_distr_cvec_eval_pdf( const double *x, struct unur_distr *distr )
500 /*----------------------------------------------------------------------*/
501 /* evaluate PDF of distribution at x */
502 /* */
503 /* parameters: */
504 /* x ... argument for pdf */
505 /* distr ... pointer to distribution object */
506 /* */
507 /* return: */
508 /* PDF(x) */
509 /*----------------------------------------------------------------------*/
510 {
511 /* check arguments */
512 _unur_check_NULL( NULL, distr, INFINITY );
513 _unur_check_distr_object( distr, CVEC, INFINITY );
514
515 if (DISTR.pdf == NULL) {
516 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
517 return INFINITY;
518 }
519
520 return _unur_cvec_PDF(x,distr);
521 } /* end of unur_distr_cvec_eval_pdf() */
522
523 /*---------------------------------------------------------------------------*/
524
525 int
unur_distr_cvec_eval_dpdf(double * result,const double * x,struct unur_distr * distr)526 unur_distr_cvec_eval_dpdf( double *result, const double *x, struct unur_distr *distr )
527 /*----------------------------------------------------------------------*/
528 /* evaluate gradient of PDF of distribution at x */
529 /* */
530 /* parameters: */
531 /* result ... to store grad (PDF(x)) */
532 /* x ... argument for dPDF */
533 /* distr ... pointer to distribution object */
534 /* */
535 /* return: */
536 /* UNUR_SUCCESS ... on success */
537 /* error code ... on error */
538 /*----------------------------------------------------------------------*/
539 {
540 /* check arguments */
541 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
542 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
543
544 if (DISTR.dpdf == NULL) {
545 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
546 return UNUR_ERR_DISTR_DATA;
547 }
548
549 return _unur_cvec_dPDF(result,x,distr);
550 } /* end of unur_distr_cvec_eval_dpdf() */
551
552 /*---------------------------------------------------------------------------*/
553
554 double
unur_distr_cvec_eval_pdpdf(const double * x,int coord,struct unur_distr * distr)555 unur_distr_cvec_eval_pdpdf( const double *x, int coord, struct unur_distr *distr )
556 /*----------------------------------------------------------------------*/
557 /* evaluate partial derivative of PDF of distribution at x for */
558 /* coordinate coord. */
559 /* */
560 /* parameters: */
561 /* x ... argument for pdf */
562 /* coord ... coordinate for partial derivative */
563 /* distr ... pointer to distribution object */
564 /* */
565 /* return: */
566 /* dPDF(x) / d(coord) */
567 /*----------------------------------------------------------------------*/
568 {
569 /* check arguments */
570 _unur_check_NULL( NULL, distr, INFINITY );
571 _unur_check_distr_object( distr, CVEC, INFINITY );
572
573 if (DISTR.pdpdf == NULL) {
574 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
575 return INFINITY;
576 }
577
578 if (coord < 0 || coord >= distr->dim) {
579 _unur_error(distr->name,UNUR_ERR_DISTR_DOMAIN,"invalid coordinate");
580 return INFINITY;
581 }
582
583 return _unur_cvec_pdPDF(x,coord,distr);
584 } /* end of unur_distr_cvec_eval_pdpdf() */
585
586 /*---------------------------------------------------------------------------*/
587
588 int
unur_distr_cvec_set_logpdf(struct unur_distr * distr,UNUR_FUNCT_CVEC * logpdf)589 unur_distr_cvec_set_logpdf( struct unur_distr *distr, UNUR_FUNCT_CVEC *logpdf )
590 /*----------------------------------------------------------------------*/
591 /* set logPDF of distribution */
592 /* */
593 /* parameters: */
594 /* distr ... pointer to distribution object */
595 /* logpdf ... pointer to logPDF */
596 /* */
597 /* return: */
598 /* UNUR_SUCCESS ... on success */
599 /* error code ... on error */
600 /*----------------------------------------------------------------------*/
601 {
602 /* check arguments */
603 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
604 _unur_check_NULL( distr->name, logpdf, UNUR_ERR_NULL);
605 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
606
607 /* we do not allow overwriting a PDF */
608 if (DISTR.pdf != NULL || DISTR.logpdf != NULL) {
609 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of logPDF not allowed");
610 return UNUR_ERR_DISTR_SET;
611 }
612
613 /* changelog */
614 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
615 /* derived parameters like mode, area, etc. might be wrong now! */
616
617 DISTR.logpdf = logpdf;
618 DISTR.pdf = _unur_distr_cvec_eval_pdf_from_logpdf;
619
620 return UNUR_SUCCESS;
621
622 } /* end of unur_distr_cvec_set_logpdf() */
623
624 /*---------------------------------------------------------------------------*/
625
626 double
_unur_distr_cvec_eval_pdf_from_logpdf(const double * x,struct unur_distr * distr)627 _unur_distr_cvec_eval_pdf_from_logpdf( const double *x, struct unur_distr *distr )
628 /*----------------------------------------------------------------------*/
629 /* evaluate PDF of distribution at x */
630 /* wrapper when only logPDF is given */
631 /* */
632 /* parameters: */
633 /* x ... argument for pdf */
634 /* distr ... pointer to distribution object */
635 /* */
636 /* return: */
637 /* PDF(x) */
638 /*----------------------------------------------------------------------*/
639 {
640 if (DISTR.logpdf == NULL) {
641 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
642 return INFINITY;
643 }
644
645 return exp(_unur_cvec_logPDF(x,distr));
646 } /* end of _unur_distr_cvec_eval_pdf_from_logpdf() */
647
648 /*---------------------------------------------------------------------------*/
649
650 int
unur_distr_cvec_set_dlogpdf(struct unur_distr * distr,UNUR_VFUNCT_CVEC * dlogpdf)651 unur_distr_cvec_set_dlogpdf( struct unur_distr *distr, UNUR_VFUNCT_CVEC *dlogpdf )
652 /*----------------------------------------------------------------------*/
653 /* set gradient of logPDF of distribution */
654 /* */
655 /* parameters: */
656 /* distr ... pointer to distribution object */
657 /* dlogpdf ... pointer to gradient of logPDF */
658 /* */
659 /* return: */
660 /* UNUR_SUCCESS ... on success */
661 /* error code ... on error */
662 /*----------------------------------------------------------------------*/
663 {
664 /* check arguments */
665 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
666 _unur_check_NULL( distr->name, dlogpdf, UNUR_ERR_NULL );
667 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
668
669 /* we do not allow overwriting a dlogPDF */
670 if (DISTR.dpdf != NULL || DISTR.dlogpdf != NULL) {
671 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of dlogPDF not allowed");
672 return UNUR_ERR_DISTR_SET;
673 }
674
675 /* changelog */
676 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
677 /* derived parameters like mode, area, etc. might be wrong now! */
678
679 DISTR.dlogpdf = dlogpdf;
680 DISTR.dpdf = _unur_distr_cvec_eval_dpdf_from_dlogpdf;
681
682 return UNUR_SUCCESS;
683 } /* end of _unur_distr_cvec_set_dlogpdf() */
684
685 /*---------------------------------------------------------------------------*/
686
687 int
_unur_distr_cvec_eval_dpdf_from_dlogpdf(double * result,const double * x,struct unur_distr * distr)688 _unur_distr_cvec_eval_dpdf_from_dlogpdf( double *result, const double *x, struct unur_distr *distr )
689 /*----------------------------------------------------------------------*/
690 /* evaluate gradient of PDF of distribution at x */
691 /* wrapper when only gradient of logPDF is given */
692 /* */
693 /* parameters: */
694 /* result ... to store grad (PDF(x)) */
695 /* x ... argument for dPDF */
696 /* distr ... pointer to distribution object */
697 /* */
698 /* return: */
699 /* UNUR_SUCCESS ... on success */
700 /* error code ... on error */
701 /*----------------------------------------------------------------------*/
702 {
703 int ret, i;
704 double fx;
705
706 if (DISTR.logpdf == NULL || DISTR.dlogpdf == NULL) {
707 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
708 return UNUR_ERR_DISTR_DATA;
709 }
710
711 fx = exp(unur_distr_cvec_eval_logpdf( x, distr ));
712 if (!_unur_isfinite(fx)) return UNUR_ERR_DISTR_DATA;
713
714 ret = _unur_cvec_dlogPDF(result,x,distr);
715 for (i=0; i<distr->dim; i++)
716 result[i] *= fx;
717
718 return ret;
719 } /* end of _unur_distr_cvec_eval_dpdf_from_dlogpdf() */
720
721 /*---------------------------------------------------------------------------*/
722
723 int
unur_distr_cvec_set_pdlogpdf(struct unur_distr * distr,UNUR_FUNCTD_CVEC * pdlogpdf)724 unur_distr_cvec_set_pdlogpdf( struct unur_distr *distr, UNUR_FUNCTD_CVEC *pdlogpdf )
725 /*----------------------------------------------------------------------*/
726 /* set partial derivative of logPDF of distribution */
727 /* */
728 /* parameters: */
729 /* distr ... pointer to distribution object */
730 /* pdlogpdf ... pointer to partial derivative of logPDF */
731 /* */
732 /* return: */
733 /* UNUR_SUCCESS ... on success */
734 /* error code ... on error */
735 /*----------------------------------------------------------------------*/
736 {
737 /* check arguments */
738 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
739 _unur_check_NULL( distr->name, pdlogpdf, UNUR_ERR_NULL);
740 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
741
742 /* we do not allow overwriting a PDF */
743 if (DISTR.pdpdf != NULL || DISTR.pdlogpdf != NULL) {
744 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of pdlogPDF not allowed");
745 return UNUR_ERR_DISTR_SET;
746 }
747
748 /* changelog */
749 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
750 /* derived parameters like mode, area, etc. might be wrong now! */
751
752 DISTR.pdlogpdf = pdlogpdf;
753 DISTR.pdpdf = _unur_distr_cvec_eval_pdpdf_from_pdlogpdf;
754
755 return UNUR_SUCCESS;
756
757 } /* end of unur_distr_cvec_set_pdlogpdf() */
758
759 /*---------------------------------------------------------------------------*/
760
761 double
_unur_distr_cvec_eval_pdpdf_from_pdlogpdf(const double * x,int coord,struct unur_distr * distr)762 _unur_distr_cvec_eval_pdpdf_from_pdlogpdf( const double *x, int coord, struct unur_distr *distr )
763 /*----------------------------------------------------------------------*/
764 /* evaluate partial derivative of PDF of distribution at x */
765 /* wrapper when only logPDF is given */
766 /* */
767 /* parameters: */
768 /* x ... argument for pdf */
769 /* distr ... pointer to distribution object */
770 /* */
771 /* return: */
772 /* PDF(x) */
773 /*----------------------------------------------------------------------*/
774 {
775 double fx;
776
777 if (DISTR.logpdf == NULL || DISTR.pdlogpdf == NULL) {
778 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
779 return INFINITY;
780 }
781
782 if (coord < 0 || coord >= distr->dim) {
783 _unur_error(distr->name,UNUR_ERR_DISTR_DOMAIN,"invalid coordinate");
784 return INFINITY;
785 }
786
787 fx = exp(unur_distr_cvec_eval_logpdf( x, distr ));
788 if (!_unur_isfinite(fx)) return INFINITY;
789
790 return fx * _unur_cvec_pdlogPDF(x,coord,distr);
791 } /* end of _unur_distr_cvec_eval_pdpdf_from_pdlogpdf() */
792
793 /*---------------------------------------------------------------------------*/
794
795 UNUR_FUNCT_CVEC *
unur_distr_cvec_get_logpdf(const struct unur_distr * distr)796 unur_distr_cvec_get_logpdf( const struct unur_distr *distr )
797 /*----------------------------------------------------------------------*/
798 /* get pointer to logPDF of distribution */
799 /* */
800 /* parameters: */
801 /* distr ... pointer to distribution object */
802 /* */
803 /* return: */
804 /* pointer to logPDF */
805 /*----------------------------------------------------------------------*/
806 {
807 /* check arguments */
808 _unur_check_NULL( NULL, distr, NULL );
809 _unur_check_distr_object( distr, CVEC, NULL );
810
811 return DISTR.logpdf;
812 } /* end of unur_distr_cvec_get_logpdf() */
813
814 /*---------------------------------------------------------------------------*/
815
816 UNUR_VFUNCT_CVEC *
unur_distr_cvec_get_dlogpdf(const struct unur_distr * distr)817 unur_distr_cvec_get_dlogpdf( const struct unur_distr *distr )
818 /*----------------------------------------------------------------------*/
819 /* get pointer to gradient of logPDF of distribution */
820 /* */
821 /* parameters: */
822 /* distr ... pointer to distribution object */
823 /* */
824 /* return: */
825 /* pointer to gradient of logPDF */
826 /*----------------------------------------------------------------------*/
827 {
828 /* check arguments */
829 _unur_check_NULL( NULL, distr, NULL );
830 _unur_check_distr_object( distr, CVEC, NULL );
831
832 return DISTR.dlogpdf;
833 } /* end of unur_distr_cvec_get_dlogpdf() */
834
835 /*---------------------------------------------------------------------------*/
836
837 UNUR_FUNCTD_CVEC *
unur_distr_cvec_get_pdlogpdf(const struct unur_distr * distr)838 unur_distr_cvec_get_pdlogpdf( const struct unur_distr *distr )
839 /*----------------------------------------------------------------------*/
840 /* get pointer to partial derivative of logPDF of distribution */
841 /* */
842 /* parameters: */
843 /* distr ... pointer to distribution object */
844 /* */
845 /* return: */
846 /* pointer to partial derivative of logPDF */
847 /*----------------------------------------------------------------------*/
848 {
849 /* check arguments */
850 _unur_check_NULL( NULL, distr, NULL );
851 _unur_check_distr_object( distr, CVEC, NULL );
852
853 return DISTR.pdlogpdf;
854 } /* end of unur_distr_cvec_get_pdlogpdf() */
855
856 /*---------------------------------------------------------------------------*/
857
858 double
unur_distr_cvec_eval_logpdf(const double * x,struct unur_distr * distr)859 unur_distr_cvec_eval_logpdf( const double *x, struct unur_distr *distr )
860 /*----------------------------------------------------------------------*/
861 /* evaluate logPDF of distribution at x */
862 /* */
863 /* parameters: */
864 /* x ... argument for logpdf */
865 /* distr ... pointer to distribution object */
866 /* */
867 /* return: */
868 /* logPDF(x) */
869 /*----------------------------------------------------------------------*/
870 {
871 /* check arguments */
872 _unur_check_NULL( NULL, distr, INFINITY );
873 _unur_check_distr_object( distr, CVEC, INFINITY );
874
875 if (DISTR.logpdf == NULL) {
876 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
877 return INFINITY;
878 }
879
880 return _unur_cvec_logPDF(x,distr);
881 } /* end of unur_distr_cvec_eval_logpdf() */
882
883 /*---------------------------------------------------------------------------*/
884
885 int
unur_distr_cvec_eval_dlogpdf(double * result,const double * x,struct unur_distr * distr)886 unur_distr_cvec_eval_dlogpdf( double *result, const double *x, struct unur_distr *distr )
887 /*----------------------------------------------------------------------*/
888 /* evaluate gradient of logPDF of distribution at x */
889 /* */
890 /* parameters: */
891 /* result ... to store grad (logPDF(x)) */
892 /* x ... argument for dlogPDF */
893 /* distr ... pointer to distribution object */
894 /* */
895 /* return: */
896 /* UNUR_SUCCESS ... on success */
897 /* error code ... on error */
898 /*----------------------------------------------------------------------*/
899 {
900 /* check arguments */
901 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
902 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
903
904 if (DISTR.dlogpdf == NULL) {
905 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
906 return UNUR_ERR_DISTR_DATA;
907 }
908
909 return _unur_cvec_dlogPDF(result,x,distr);
910 } /* end of unur_distr_cvec_eval_dlogpdf() */
911
912 /*---------------------------------------------------------------------------*/
913
914 double
unur_distr_cvec_eval_pdlogpdf(const double * x,int coord,struct unur_distr * distr)915 unur_distr_cvec_eval_pdlogpdf( const double *x, int coord, struct unur_distr *distr )
916 /*----------------------------------------------------------------------*/
917 /* evaluate of partial derivative of logPDF of distribution at x */
918 /* in direction coord */
919 /* */
920 /* parameters: */
921 /* x ... argument for logpdf */
922 /* coord ... coordinate for partial derivative */
923 /* distr ... pointer to distribution object */
924 /* */
925 /* return: */
926 /* d logPDF(x) / d(coord) */
927 /*----------------------------------------------------------------------*/
928 {
929 /* check arguments */
930 _unur_check_NULL( NULL, distr, INFINITY );
931 _unur_check_distr_object( distr, CVEC, INFINITY );
932
933 if (DISTR.pdlogpdf == NULL) {
934 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
935 return INFINITY;
936 }
937
938 if (coord < 0 || coord >= distr->dim) {
939 _unur_error(distr->name,UNUR_ERR_DISTR_DOMAIN,"invalid coordinate");
940 return INFINITY;
941 }
942
943 return _unur_cvec_pdlogPDF(x,coord,distr);
944 } /* end of unur_distr_cvec_eval_pdlogpdf() */
945
946 /*---------------------------------------------------------------------------*/
947
948 int
unur_distr_cvec_set_domain_rect(struct unur_distr * distr,const double * lowerleft,const double * upperright)949 unur_distr_cvec_set_domain_rect( struct unur_distr *distr, const double *lowerleft, const double *upperright )
950 /*----------------------------------------------------------------------*/
951 /* Set rectangular domain */
952 /* */
953 /* parameters: */
954 /* distr ... pointer to distribution object */
955 /* lowerleft ... lower left vertex of rectanlge */
956 /* upperright ... upper right vertex of rectanlge */
957 /* */
958 /* return: */
959 /* UNUR_SUCCESS ... on success */
960 /* error code ... on error */
961 /*----------------------------------------------------------------------*/
962 {
963 int i;
964
965 /* check arguments */
966 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
967 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
968 _unur_check_NULL( distr->name, lowerleft, UNUR_ERR_NULL );
969 _unur_check_NULL( distr->name, upperright, UNUR_ERR_NULL );
970
971 /* check new parameter for distribution */
972 for (i=0; i<distr->dim; i++) {
973 if (!(lowerleft[i] < upperright[i] * (1.-UNUR_SQRT_DBL_EPSILON))) {
974 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"domain, left >= right");
975 return UNUR_ERR_DISTR_SET;
976 }
977 }
978
979 /* store data */
980 DISTR.domainrect = _unur_xrealloc(DISTR.domainrect, 2 * distr->dim * sizeof(double));
981 for (i=0; i<distr->dim; i++) {
982 DISTR.domainrect[2*i] = lowerleft[i];
983 DISTR.domainrect[2*i+1] = upperright[i];
984 }
985
986 /* changelog */
987 distr->set |= UNUR_DISTR_SET_DOMAIN | UNUR_DISTR_SET_DOMAINBOUNDED;
988 /* we silently assume here that at least one of the given coordinates */
989 /* is finite. */
990
991 /* we have to mark all derived parameters as unknown */
992 distr->set &= ~(UNUR_DISTR_SET_STDDOMAIN |
993 UNUR_DISTR_SET_MASK_DERIVED );
994
995 if (distr->base) {
996 /* for derived distributions (e.g. order statistics)
997 we also have to mark derived parameters as unknown */
998 distr->base->set &= ~(UNUR_DISTR_SET_STDDOMAIN |
999 UNUR_DISTR_SET_MASK_DERIVED );
1000 /* We only can set the domain of the base distribution if it */
1001 /* is of the same type as distr */
1002 if ( distr->base->type == UNUR_DISTR_CVEC ) {
1003 if (unur_distr_cvec_set_domain_rect(distr->base, lowerleft, upperright)!=UNUR_SUCCESS)
1004 return UNUR_ERR_DISTR_SET;
1005 }
1006 }
1007
1008 /* o.k. */
1009 return UNUR_SUCCESS;
1010 } /* end unur_distr_cvec_set_domain_rect() */
1011
1012 /*---------------------------------------------------------------------------*/
1013
1014 int
_unur_distr_cvec_has_boundeddomain(const struct unur_distr * distr)1015 _unur_distr_cvec_has_boundeddomain( const struct unur_distr *distr )
1016 /*----------------------------------------------------------------------*/
1017 /* Check whether distr has a bounded domain */
1018 /* */
1019 /* parameters: */
1020 /* distr ... pointer to distribution object */
1021 /* */
1022 /* return: */
1023 /* TRUE ... if domain is bounded */
1024 /* FALSE ... otherwise or in case of an error */
1025 /*----------------------------------------------------------------------*/
1026 {
1027 int i;
1028 double *domain;
1029
1030 /* check arguments */
1031 CHECK_NULL( distr, FALSE );
1032 COOKIE_CHECK(distr,CK_DISTR_CVEC,FALSE);
1033
1034 if (! (distr->set & UNUR_DISTR_SET_DOMAINBOUNDED &&
1035 DISTR.domainrect))
1036 return FALSE;
1037
1038 domain = DISTR.domainrect;
1039 for (i=0; i < 2*distr->dim; i++)
1040 if (!_unur_isfinite(domain[i]))
1041 return FALSE;
1042
1043 return TRUE;
1044 } /* end of _unur_distr_cvec_has_boundeddomain() */
1045
1046 /*---------------------------------------------------------------------------*/
1047
1048 int
_unur_distr_cvec_is_indomain(const double * x,const struct unur_distr * distr)1049 _unur_distr_cvec_is_indomain( const double *x, const struct unur_distr *distr )
1050 /*----------------------------------------------------------------------*/
1051 /* Check whether x falls into domain */
1052 /* */
1053 /* parameters: */
1054 /* x ... pointer to point */
1055 /* distr ... pointer to distribution object */
1056 /* */
1057 /* return: */
1058 /* TRUE ... if x in domain */
1059 /* FALSE ... otherwise or in case of an error */
1060 /*----------------------------------------------------------------------*/
1061 {
1062 int i;
1063 double *domain;
1064
1065 /* check arguments */
1066 CHECK_NULL( distr, FALSE );
1067 COOKIE_CHECK(distr,CK_DISTR_CVEC,FALSE);
1068
1069 domain = DISTR.domainrect;
1070 if (domain==NULL)
1071 /* unbounded domain */
1072 return TRUE;
1073
1074 for (i=0; i<distr->dim; i++) {
1075 if (x[i] < domain[2*i] || x[i] > domain[2*i+1])
1076 return FALSE;
1077 }
1078
1079 return TRUE;
1080 } /* end of _unur_distr_cvec_is_indomain() */
1081
1082 int
unur_distr_cvec_is_indomain(const double * x,const struct unur_distr * distr)1083 unur_distr_cvec_is_indomain( const double *x, const struct unur_distr *distr )
1084 {
1085 /* check arguments */
1086 _unur_check_NULL( NULL, distr, FALSE );
1087 _unur_check_distr_object( distr, CVEC, FALSE );
1088
1089 return _unur_distr_cvec_is_indomain(x, distr);
1090 } /* end of unur_distr_cvec_is_indomain() */
1091
1092 /*---------------------------------------------------------------------------*/
1093
1094 int
unur_distr_cvec_set_mean(struct unur_distr * distr,const double * mean)1095 unur_distr_cvec_set_mean( struct unur_distr *distr, const double *mean )
1096 /*----------------------------------------------------------------------*/
1097 /* set mean vector of distribution */
1098 /* */
1099 /* parameters: */
1100 /* distr ... pointer to distribution object */
1101 /* mean ... mean vector of distribution */
1102 /* */
1103 /* return: */
1104 /* UNUR_SUCCESS ... on success */
1105 /* error code ... on error */
1106 /*----------------------------------------------------------------------*/
1107 {
1108 int i;
1109
1110 /* check arguments */
1111 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1112 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1113
1114 /* we have to allocate memory first */
1115 if (DISTR.mean == NULL)
1116 DISTR.mean = _unur_xmalloc( distr->dim * sizeof(double) );
1117
1118 if (mean)
1119 /* mean vector given --> copy */
1120 memcpy( DISTR.mean, mean, distr->dim * sizeof(double) );
1121
1122 else /* mean == NULL --> use zero vector instead */
1123 for (i=0; i<distr->dim; i++)
1124 DISTR.mean[i] = 0.;
1125
1126 /* changelog */
1127 distr->set |= UNUR_DISTR_SET_MEAN;
1128
1129 /* o.k. */
1130 return UNUR_SUCCESS;
1131 } /* end of unur_distr_cvec_set_mean() */
1132
1133 /*---------------------------------------------------------------------------*/
1134
1135 const double *
unur_distr_cvec_get_mean(const struct unur_distr * distr)1136 unur_distr_cvec_get_mean( const struct unur_distr *distr )
1137 /*----------------------------------------------------------------------*/
1138 /* get mean vector of distribution */
1139 /* */
1140 /* parameters: */
1141 /* distr ... pointer to distribution object */
1142 /* */
1143 /* return: */
1144 /* pointer to mean of distribution */
1145 /* */
1146 /* error: */
1147 /* return NULL */
1148 /*----------------------------------------------------------------------*/
1149 {
1150 /* check arguments */
1151 _unur_check_NULL( NULL, distr, NULL );
1152 _unur_check_distr_object( distr, CVEC, NULL );
1153
1154 /* mean vector known ? */
1155 if ( !(distr->set & UNUR_DISTR_SET_MEAN) ) {
1156 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"mean");
1157 return NULL;
1158 }
1159
1160 return DISTR.mean;
1161
1162 } /* end of unur_distr_cvec_get_mean() */
1163
1164 /*---------------------------------------------------------------------------*/
1165
1166 int
unur_distr_cvec_set_covar(struct unur_distr * distr,const double * covar)1167 unur_distr_cvec_set_covar( struct unur_distr *distr, const double *covar )
1168 /*----------------------------------------------------------------------*/
1169 /* set covariance matrix of distribution. */
1170 /* as a side effect it also computes its cholesky factor. */
1171 /* */
1172 /* parameters: */
1173 /* distr ... pointer to distribution object */
1174 /* covar ... covariance matrix of distribution */
1175 /* */
1176 /* return: */
1177 /* UNUR_SUCCESS ... on success */
1178 /* error code ... on error */
1179 /*----------------------------------------------------------------------*/
1180 {
1181 #define idx(a,b) ((a)*dim+(b))
1182
1183 int i,j;
1184 int dim;
1185
1186 /* check arguments */
1187 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1188 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1189
1190 dim = distr->dim;
1191
1192 /* mark as unknown */
1193 distr->set &= ~( UNUR_DISTR_SET_COVAR
1194 | UNUR_DISTR_SET_COVAR_IDENT
1195 | UNUR_DISTR_SET_CHOLESKY
1196 | UNUR_DISTR_SET_COVAR_INV );
1197
1198 /* we have to allocate memory first */
1199 if (DISTR.covar == NULL)
1200 DISTR.covar = _unur_xmalloc( dim * dim * sizeof(double) );
1201 if (DISTR.cholesky == NULL)
1202 DISTR.cholesky = _unur_xmalloc( dim * dim * sizeof(double) );
1203
1204 /* if covar == NULL --> use identity matrix */
1205 if (covar==NULL) {
1206 for (i=0; i<dim; i++) {
1207 for (j=0; j<dim; j++) {
1208 DISTR.covar[idx(i,j)] = (i==j) ? 1. : 0.;
1209 DISTR.cholesky[idx(i,j)] = (i==j) ? 1. : 0.;
1210 }
1211 }
1212 distr->set |= UNUR_DISTR_SET_COVAR_IDENT;
1213 }
1214
1215 /* covariance matrix given --> copy data */
1216 else {
1217
1218 /* check covariance matrix: diagonal entries > 0 */
1219 for (i=0; i<dim*dim; i+= dim+1)
1220 if (covar[i] <= 0.) {
1221 _unur_error(distr->name ,UNUR_ERR_DISTR_DOMAIN,"variance <= 0");
1222 return UNUR_ERR_DISTR_DOMAIN;
1223 }
1224
1225 /* check for symmetry */
1226 for (i=0; i<dim; i++)
1227 for (j=i+1; j<dim; j++)
1228 if (!_unur_FP_same(covar[i*dim+j],covar[j*dim+i])) {
1229 _unur_error(distr->name ,UNUR_ERR_DISTR_DOMAIN,
1230 "covariance matrix not symmetric");
1231 return UNUR_ERR_DISTR_DOMAIN;
1232 }
1233
1234 /* copy data */
1235 memcpy( DISTR.covar, covar, dim * dim * sizeof(double) );
1236
1237 /* compute Cholesky decomposition and check for positive definitness */
1238 if (_unur_matrix_cholesky_decomposition(dim, covar, DISTR.cholesky) != UNUR_SUCCESS) {
1239 _unur_error(distr->name, UNUR_ERR_DISTR_DOMAIN,
1240 "covariance matrix not positive definite");
1241 return UNUR_ERR_DISTR_DOMAIN;
1242 }
1243
1244 }
1245
1246 /* changelog */
1247 distr->set |= UNUR_DISTR_SET_COVAR | UNUR_DISTR_SET_CHOLESKY;
1248
1249 /* o.k. */
1250 return UNUR_SUCCESS;
1251
1252 #undef idx
1253 } /* end of unur_distr_cvec_set_covar() */
1254
1255 /*---------------------------------------------------------------------------*/
1256
1257 int
unur_distr_cvec_set_covar_inv(struct unur_distr * distr,const double * covar_inv)1258 unur_distr_cvec_set_covar_inv( struct unur_distr *distr, const double *covar_inv )
1259 /*----------------------------------------------------------------------*/
1260 /* set inverse of covariance matrix of distribution. */
1261 /* */
1262 /* parameters: */
1263 /* distr ... pointer to distribution object */
1264 /* covar_inv ... inverse of covariance matrix of distribution */
1265 /* */
1266 /* return: */
1267 /* UNUR_SUCCESS ... on success */
1268 /* error code ... on error */
1269 /*----------------------------------------------------------------------*/
1270 {
1271 #define idx(a,b) ((a)*dim+(b))
1272
1273 int i,j;
1274 int dim;
1275
1276 /* check arguments */
1277 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1278 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1279
1280 dim = distr->dim;
1281
1282 /* mark as unknown */
1283 distr->set &= ~(UNUR_DISTR_SET_COVAR_INV);
1284
1285 /* we have to allocate memory first */
1286 if (DISTR.covar_inv == NULL)
1287 DISTR.covar_inv = _unur_xmalloc( dim * dim * sizeof(double) );
1288
1289 /* if covar_inv == NULL --> use identity matrix */
1290 if (covar_inv==NULL)
1291 for (i=0; i<dim; i++)
1292 for (j=0; j<dim; j++)
1293 DISTR.covar_inv[idx(i,j)] = (i==j) ? 1. : 0.;
1294
1295 /* inverse of covariance matrix given --> copy data */
1296 else {
1297
1298 /* check inverse of covariance matrix: diagonal entries > 0 */
1299 for (i=0; i<dim*dim; i+= dim+1)
1300 if (covar_inv[i] <= 0.) {
1301 _unur_error(distr->name ,UNUR_ERR_DISTR_DOMAIN,"diagonals <= 0");
1302 return UNUR_ERR_DISTR_DOMAIN;
1303 }
1304
1305 /* check for symmetry */
1306 for (i=0; i<dim; i++)
1307 for (j=i+1; j<dim; j++)
1308 if (!_unur_FP_same(covar_inv[i*dim+j],covar_inv[j*dim+i])) {
1309 _unur_error(distr->name ,UNUR_ERR_DISTR_DOMAIN,
1310 "inverse of covariance matrix not symmetric");
1311 return UNUR_ERR_DISTR_DOMAIN;
1312 }
1313
1314 /* copy data */
1315 memcpy( DISTR.covar_inv, covar_inv, dim * dim * sizeof(double) );
1316
1317 }
1318
1319 /* changelog */
1320 distr->set |= UNUR_DISTR_SET_COVAR_INV;
1321
1322 /* o.k. */
1323 return UNUR_SUCCESS;
1324
1325 #undef idx
1326 } /* end of unur_distr_cvec_set_covar_inv() */
1327
1328 /*---------------------------------------------------------------------------*/
1329
1330 const double *
unur_distr_cvec_get_covar(const struct unur_distr * distr)1331 unur_distr_cvec_get_covar( const struct unur_distr *distr )
1332 /*----------------------------------------------------------------------*/
1333 /* get covariance matrix of distribution */
1334 /* */
1335 /* parameters: */
1336 /* distr ... pointer to distribution object */
1337 /* */
1338 /* return: */
1339 /* pointer to covariance matrix of distribution */
1340 /* */
1341 /* error: */
1342 /* return NULL */
1343 /*----------------------------------------------------------------------*/
1344 {
1345 /* check arguments */
1346 _unur_check_NULL( NULL, distr, NULL );
1347 _unur_check_distr_object( distr, CVEC, NULL );
1348
1349 /* covariance matrix known ? */
1350 if ( !(distr->set & UNUR_DISTR_SET_COVAR) ) {
1351 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"covariance matrix");
1352 return NULL;
1353 }
1354
1355 return DISTR.covar;
1356
1357 } /* end of unur_distr_cvec_get_covar() */
1358
1359 /*---------------------------------------------------------------------------*/
1360
1361 const double *
unur_distr_cvec_get_cholesky(const struct unur_distr * distr)1362 unur_distr_cvec_get_cholesky( const struct unur_distr *distr )
1363 /*----------------------------------------------------------------------*/
1364 /* get cholesky factor of the covariance matrix of distribution */
1365 /* */
1366 /* parameters: */
1367 /* distr ... pointer to distribution object */
1368 /* */
1369 /* return: */
1370 /* pointer to cholesky factor */
1371 /*----------------------------------------------------------------------*/
1372 {
1373 /* check arguments */
1374 _unur_check_NULL( NULL, distr, NULL );
1375 _unur_check_distr_object( distr, CVEC, NULL );
1376
1377 /* covariance matrix known ? */
1378 if ( !(distr->set & UNUR_DISTR_SET_CHOLESKY) ) {
1379 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"covariance matrix");
1380 return NULL;
1381 }
1382
1383 return DISTR.cholesky;
1384
1385 } /* end of unur_distr_cvec_get_covar_cholesky() */
1386
1387 /*---------------------------------------------------------------------------*/
1388
1389 const double *
unur_distr_cvec_get_covar_inv(struct unur_distr * distr)1390 unur_distr_cvec_get_covar_inv ( struct unur_distr *distr )
1391 /*----------------------------------------------------------------------*/
1392 /* get inverse covariance matrix of distribution */
1393 /* */
1394 /* parameters: */
1395 /* distr ... pointer to distribution object */
1396 /* */
1397 /* return: */
1398 /* pointer to inverse of covariance matrix */
1399 /*----------------------------------------------------------------------*/
1400 {
1401
1402 double det; /* determinant of covariance matrix */
1403 int dim;
1404
1405 /* check arguments */
1406 _unur_check_NULL( NULL, distr, NULL );
1407 _unur_check_distr_object( distr, CVEC, NULL );
1408
1409 dim = distr->dim;
1410
1411 /* covariance matrix known ? */
1412 if ( !(distr->set & UNUR_DISTR_SET_COVAR) ) {
1413 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"covariance matrix not known");
1414 return NULL;
1415 }
1416
1417 /* allocate memory */
1418 if (DISTR.covar_inv == NULL)
1419 DISTR.covar_inv = _unur_xmalloc( dim * dim * sizeof(double) );
1420
1421 if ( !(distr->set & UNUR_DISTR_SET_COVAR_INV) ) {
1422 /* calculate inverse covariance matrix */
1423 if (_unur_matrix_invert_matrix(dim, DISTR.covar, DISTR.covar_inv, &det) != UNUR_SUCCESS) {
1424 _unur_error(distr->name ,UNUR_ERR_DISTR_DOMAIN,"cannot compute inverse of covariance");
1425 return NULL;
1426 }
1427 }
1428
1429 /* changelog */
1430 distr->set |= UNUR_DISTR_SET_COVAR_INV;
1431
1432 return DISTR.covar_inv;
1433
1434 } /* end of unur_distr_cvec_get_covar_inv() */
1435
1436 /*---------------------------------------------------------------------------*/
1437
1438 int
unur_distr_cvec_set_rankcorr(struct unur_distr * distr,const double * rankcorr)1439 unur_distr_cvec_set_rankcorr( struct unur_distr *distr, const double *rankcorr )
1440 /*----------------------------------------------------------------------*/
1441 /* Set rank-correlation matrix of distribution. */
1442 /* The given matrix is checked for symmetry and positive definitness. */
1443 /* The diagonal entries must be equal to 1. */
1444 /* As a side effect it also computes its cholesky factor. */
1445 /* */
1446 /* parameters: */
1447 /* distr ... pointer to distribution object */
1448 /* rankcorr ... rankcorrelation matrix of distribution */
1449 /* */
1450 /* return: */
1451 /* UNUR_SUCCESS ... on success */
1452 /* error code ... on error */
1453 /*----------------------------------------------------------------------*/
1454 {
1455 #define idx(a,b) ((a)*dim+(b))
1456
1457 int i,j;
1458 int dim;
1459
1460 /* check arguments */
1461 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1462 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1463
1464 dim = distr->dim;
1465
1466 /* mark as unknown */
1467 distr->set &= ~(UNUR_DISTR_SET_RANKCORR | UNUR_DISTR_SET_RK_CHOLESKY);
1468
1469 /* we have to allocate memory first */
1470 if (DISTR.rankcorr == NULL)
1471 DISTR.rankcorr = _unur_xmalloc( dim * dim * sizeof(double) );
1472 if (DISTR.rk_cholesky == NULL)
1473 DISTR.rk_cholesky = _unur_xmalloc( dim * dim * sizeof(double) );
1474
1475 /* if rankcorr == NULL --> use identity matrix */
1476 if (rankcorr==NULL) {
1477 for (i=0; i<dim; i++)
1478 for (j=0; j<dim; j++) {
1479 DISTR.rankcorr[idx(i,j)] = (i==j) ? 1. : 0.;
1480 DISTR.rk_cholesky[idx(i,j)] = (i==j) ? 1. : 0.;
1481 }
1482 }
1483
1484 /* rankcorriance matrix given --> copy data */
1485 else {
1486
1487 /* check rankcorriance matrix: diagonal entries == 1 */
1488 for (i=0; i<dim*dim; i+= dim+1) {
1489 if (!_unur_FP_same(rankcorr[i],1.)) {
1490 _unur_error(distr->name ,UNUR_ERR_DISTR_DOMAIN,"diagonals != 1");
1491 return UNUR_ERR_DISTR_DOMAIN;
1492 }
1493 }
1494
1495 /* check for symmetry */
1496 for (i=0; i<dim; i++)
1497 for (j=i+1; j<dim; j++)
1498 if (!_unur_FP_same(rankcorr[i*dim+j],rankcorr[j*dim+i])) {
1499 _unur_error(distr->name ,UNUR_ERR_DISTR_DOMAIN,
1500 "rank-correlation matrix not symmetric");
1501 return UNUR_ERR_DISTR_DOMAIN;
1502 }
1503
1504 /* copy data */
1505 memcpy( DISTR.rankcorr, rankcorr, dim * dim * sizeof(double) );
1506
1507 /* compute Cholesky decomposition and check for positive definitness */
1508 if (_unur_matrix_cholesky_decomposition(dim, rankcorr, DISTR.rk_cholesky) != UNUR_SUCCESS) {
1509 _unur_error(distr->name, UNUR_ERR_DISTR_DOMAIN,
1510 "rankcorriance matrix not positive definite");
1511 return UNUR_ERR_DISTR_DOMAIN;
1512 }
1513
1514 }
1515
1516 /* changelog */
1517 distr->set |= UNUR_DISTR_SET_RANKCORR | UNUR_DISTR_SET_RK_CHOLESKY;
1518
1519 /* o.k. */
1520 return UNUR_SUCCESS;
1521
1522 #undef idx
1523 } /* end of unur_distr_cvec_set_rankcorr() */
1524
1525 /*---------------------------------------------------------------------------*/
1526
1527 const double *
unur_distr_cvec_get_rankcorr(const struct unur_distr * distr)1528 unur_distr_cvec_get_rankcorr( const struct unur_distr *distr )
1529 /*----------------------------------------------------------------------*/
1530 /* get rank-correlation matrix of distribution */
1531 /* */
1532 /* parameters: */
1533 /* distr ... pointer to distribution object */
1534 /* */
1535 /* return: */
1536 /* pointer to rank-correlation matrix of distribution */
1537 /* */
1538 /* error: */
1539 /* return NULL */
1540 /*----------------------------------------------------------------------*/
1541 {
1542 /* check arguments */
1543 _unur_check_NULL( NULL, distr, NULL );
1544 _unur_check_distr_object( distr, CVEC, NULL );
1545
1546 /* rankcorriance matrix known ? */
1547 if ( !(distr->set & UNUR_DISTR_SET_RANKCORR) ) {
1548 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"rank-correlation matrix");
1549 return NULL;
1550 }
1551
1552 return DISTR.rankcorr;
1553
1554 } /* end of unur_distr_cvec_get_rankcorr() */
1555
1556 /*---------------------------------------------------------------------------*/
1557
1558 const double *
unur_distr_cvec_get_rk_cholesky(const struct unur_distr * distr)1559 unur_distr_cvec_get_rk_cholesky( const struct unur_distr *distr )
1560 /*----------------------------------------------------------------------*/
1561 /* get cholesky factor of the rank correlation matrix of distribution */
1562 /* */
1563 /* parameters: */
1564 /* distr ... pointer to distribution object */
1565 /* */
1566 /* return: */
1567 /* pointer to cholesky factor */
1568 /*----------------------------------------------------------------------*/
1569 {
1570 /* check arguments */
1571 _unur_check_NULL( NULL, distr, NULL );
1572 _unur_check_distr_object( distr, CVEC, NULL );
1573
1574 /* covariance matrix known ? */
1575 if ( !(distr->set & UNUR_DISTR_SET_RK_CHOLESKY) ) {
1576 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"rank correlation matrix");
1577 return NULL;
1578 }
1579
1580 return DISTR.rk_cholesky;
1581
1582 } /* end of unur_distr_cvec_get_rk_cholesky() */
1583
1584 /*---------------------------------------------------------------------------*/
1585
1586 int
unur_distr_cvec_set_marginals(struct unur_distr * distr,struct unur_distr * marginal)1587 unur_distr_cvec_set_marginals ( struct unur_distr *distr, struct unur_distr *marginal)
1588 /*----------------------------------------------------------------------*/
1589 /* Copy marginal distribution into distribution object. */
1590 /* Only one local copy is made. */
1591 /* */
1592 /* parameters: */
1593 /* distr ... pointer to distribution object */
1594 /* marginal ... pointer to marginal distribution object */
1595 /* */
1596 /* */
1597 /* return: */
1598 /* UNUR_SUCCESS ... on success */
1599 /* error code ... on error */
1600 /*----------------------------------------------------------------------*/
1601 {
1602 struct unur_distr *clone;
1603 int i;
1604
1605 /* check arguments */
1606 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1607 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1608 _unur_check_NULL( distr->name, marginal, UNUR_ERR_NULL );
1609 _unur_check_distr_object( marginal, CONT, UNUR_ERR_DISTR_INVALID );
1610
1611 /* first we have to check whether there is already a list of marginal distributions */
1612 if (DISTR.marginals)
1613 _unur_distr_cvec_marginals_free(DISTR.marginals, distr->dim);
1614
1615 /* make copy of marginal distribution object */
1616 clone = _unur_distr_clone( marginal );
1617
1618 /* allocate memory for array */
1619 DISTR.marginals = _unur_xmalloc (distr->dim * sizeof(struct unur_distr *));
1620
1621 /* copy pointer */
1622 for (i=0; i<distr->dim; i++)
1623 DISTR.marginals[i] = clone;
1624
1625 /* changelog */
1626 distr->set |= UNUR_DISTR_SET_MARGINAL;
1627
1628 return UNUR_SUCCESS;
1629 } /* end of unur_distr_cvec_set_marginals() */
1630
1631 /*---------------------------------------------------------------------------*/
1632
1633 int
unur_distr_cvec_set_marginal_array(struct unur_distr * distr,struct unur_distr ** marginals)1634 unur_distr_cvec_set_marginal_array ( struct unur_distr *distr, struct unur_distr **marginals)
1635 /*----------------------------------------------------------------------*/
1636 /* Copy marginal distributions into distribution object. */
1637 /* For each dimension a new copy is made even if the pointer in */
1638 /* the array marginals coincide. */
1639 /* */
1640 /* parameters: */
1641 /* distr ... pointer to distribution object */
1642 /* marginals ... pointer to array of marginal distribution objects */
1643 /* */
1644 /* */
1645 /* return: */
1646 /* UNUR_SUCCESS ... on success */
1647 /* error code ... on error */
1648 /*----------------------------------------------------------------------*/
1649 {
1650 int i;
1651
1652 /* check arguments */
1653 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1654 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1655 _unur_check_NULL( distr->name, marginals, UNUR_ERR_NULL );
1656
1657 for (i=0; i<distr->dim; i++) {
1658 _unur_check_NULL( distr->name, *(marginals+i), UNUR_ERR_NULL );
1659 _unur_check_distr_object( *(marginals+i), CONT, UNUR_ERR_DISTR_INVALID );
1660 }
1661
1662 /* first we have to check whether there is already a list of marginal distributions */
1663 if (DISTR.marginals)
1664 _unur_distr_cvec_marginals_free(DISTR.marginals, distr->dim);
1665
1666 /* allocate memory for array */
1667 DISTR.marginals = _unur_xmalloc (distr->dim * sizeof(struct unur_distr *));
1668
1669 /* make copy of marginal distribution objects */
1670 for (i=0; i<distr->dim; i++)
1671 DISTR.marginals[i] = _unur_distr_clone( *(marginals+i) );
1672
1673 /* changelog */
1674 distr->set |= UNUR_DISTR_SET_MARGINAL;
1675
1676 return UNUR_SUCCESS;
1677 } /* end of unur_distr_cvec_set_marginal_array() */
1678
1679 /*---------------------------------------------------------------------------*/
1680
1681 int
unur_distr_cvec_set_marginal_list(struct unur_distr * distr,...)1682 unur_distr_cvec_set_marginal_list ( struct unur_distr *distr, ... )
1683 /*----------------------------------------------------------------------*/
1684 /* Copy marginal distributions into distribution object. */
1685 /* For each dimenision there must be a pointer to a discribution object.*/
1686 /* */
1687 /* parameters: */
1688 /* distr ... pointer to distribution object */
1689 /* ... ... pointer to array of marginal distribution objects */
1690 /* */
1691 /* */
1692 /* return: */
1693 /* UNUR_SUCCESS ... on success */
1694 /* error code ... on error */
1695 /* */
1696 /* IMPORTANT: */
1697 /* All marginal distribution objects are destroyed after they have */
1698 /* been copied into the distribution object. */
1699 /*----------------------------------------------------------------------*/
1700 {
1701 int i;
1702 int failed = FALSE;
1703 struct unur_distr *marginal;
1704 struct unur_distr **marginal_list;
1705 va_list vargs;
1706
1707 /* check arguments */
1708 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1709 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1710
1711 /* allocate memory for array */
1712 marginal_list = _unur_xmalloc (distr->dim * sizeof(struct unur_distr *));
1713 for (i=0; i<distr->dim; i++) marginal_list[i] = NULL;
1714
1715 /* make copy of marginal distribution objects */
1716 va_start(vargs, distr);
1717 for (i=0; i<distr->dim; i++) {
1718 marginal = (struct unur_distr *) va_arg(vargs, struct unur_distr *);
1719 if (marginal) {
1720 marginal_list[i] = _unur_distr_clone( marginal );
1721 _unur_distr_free(marginal);
1722 }
1723 else {
1724 failed = TRUE;
1725 }
1726 }
1727 va_end(vargs);
1728
1729 if (failed) {
1730 /* some of the pointers are NULL pointers */
1731 _unur_distr_cvec_marginals_free(marginal_list, distr->dim);
1732 _unur_error(distr->name ,UNUR_ERR_DISTR_SET,"marginals == NULL");
1733 return UNUR_ERR_DISTR_SET;
1734 }
1735
1736 /* copy list of marginal distributions. However, first we have to check */
1737 /* whether there is already a list of marginal distributions. */
1738 if (DISTR.marginals)
1739 _unur_distr_cvec_marginals_free(DISTR.marginals, distr->dim);
1740
1741 DISTR.marginals = marginal_list;
1742
1743 /* changelog */
1744 distr->set |= UNUR_DISTR_SET_MARGINAL;
1745
1746 return UNUR_SUCCESS;
1747 } /* end of unur_distr_cvec_set_marginal_list() */
1748
1749 /*---------------------------------------------------------------------------*/
1750
1751 const struct unur_distr *
unur_distr_cvec_get_marginal(const struct unur_distr * distr,int n)1752 unur_distr_cvec_get_marginal( const struct unur_distr *distr, int n )
1753 /*----------------------------------------------------------------------*/
1754 /* Get pointer to marginal distribution object. */
1755 /* */
1756 /* parameters: */
1757 /* distr ... pointer to distribution object */
1758 /* n ... position of marginal distribution object */
1759 /* */
1760 /* return: */
1761 /* UNUR_SUCCESS ... on success */
1762 /* error code ... on error */
1763 /*----------------------------------------------------------------------*/
1764 {
1765 /* check arguments */
1766 _unur_check_NULL( NULL, distr, NULL );
1767 _unur_check_distr_object( distr, CVEC, NULL );
1768
1769 if (n<=0 || n > distr->dim) {
1770 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"n not in 1 .. dim");
1771 return NULL;
1772 }
1773
1774 /* marginal distributions known ? */
1775 if ( !(distr->set & UNUR_DISTR_SET_MARGINAL) ) {
1776 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"marginals");
1777 return NULL;
1778 }
1779
1780 _unur_check_NULL( distr->name, DISTR.marginals, NULL );
1781
1782 /* return marginal distribution object */
1783 return (DISTR.marginals[n-1]);
1784 } /* end of unur_distr_cvec_get_marginal() */
1785
1786 /*---------------------------------------------------------------------------*/
1787
1788 struct unur_distr **
_unur_distr_cvec_marginals_clone(struct unur_distr ** marginals,int dim)1789 _unur_distr_cvec_marginals_clone ( struct unur_distr **marginals, int dim )
1790 /*----------------------------------------------------------------------*/
1791 /* copy (clone) list of marginal distribution objects */
1792 /* */
1793 /* parameters: */
1794 /* marginals ... pointer to list of marginal distribution objects */
1795 /* dim ... number of marginal distributions */
1796 /* */
1797 /* return: */
1798 /* pointer to clone of list of marginal distribution objects */
1799 /*----------------------------------------------------------------------*/
1800 {
1801 struct unur_distr **clone;
1802 int i;
1803
1804 _unur_check_NULL( NULL, marginals, NULL );
1805
1806 /* check dimension */
1807 if (dim < 1) {
1808 _unur_error(NULL,UNUR_ERR_DISTR_SET,"dimension < 1");
1809 return NULL;
1810 }
1811
1812 /* allocate memory for array */
1813 clone = _unur_xmalloc (dim * sizeof(struct unur_distr *));
1814
1815 if (_unur_distr_cvec_marginals_are_equal(marginals, dim)) {
1816 clone[0] = _unur_distr_clone( marginals[0] );
1817 for (i=1; i<dim; i++)
1818 clone[i] = clone[0];
1819 }
1820
1821 else {
1822 for (i=0; i<dim; i++)
1823 clone[i] = _unur_distr_clone( marginals[i] );
1824 }
1825
1826 return clone;
1827 } /* end of _unur_distr_cvec_marginals_clone() */
1828
1829 /*---------------------------------------------------------------------------*/
1830
1831 void
_unur_distr_cvec_marginals_free(struct unur_distr ** marginals,int dim)1832 _unur_distr_cvec_marginals_free ( struct unur_distr **marginals, int dim )
1833 /*----------------------------------------------------------------------*/
1834 /* free list of marginal distribution objects */
1835 /* */
1836 /* parameters: */
1837 /* marginals ... pointer to list of marginal distribution objects */
1838 /* dim ... number of marginal distributions */
1839 /* */
1840 /* return: */
1841 /* pointer to clone of list of marginal distribution objects */
1842 /*----------------------------------------------------------------------*/
1843 {
1844 int i;
1845
1846 if (_unur_distr_cvec_marginals_are_equal(marginals,dim)) {
1847 _unur_distr_free(marginals[0]);
1848 }
1849
1850 else {
1851 for (i=0; i<dim; i++)
1852 if (marginals[i]) _unur_distr_free(marginals[i]);
1853 }
1854
1855 free (marginals);
1856 } /* end of _unur_distr_cvec_marginals_free() */
1857
1858 /*---------------------------------------------------------------------------*/
1859
1860 int
_unur_distr_cvec_marginals_are_equal(struct unur_distr ** marginals,int dim)1861 _unur_distr_cvec_marginals_are_equal( struct unur_distr **marginals, int dim )
1862 /*----------------------------------------------------------------------*/
1863 /* test whether all marginals are equal or not. */
1864 /* */
1865 /* parameters: */
1866 /* marginals ... pointer to list of marginal distribution objects */
1867 /* dim ... dimension of distribution (= number of marginals) */
1868 /* */
1869 /* return: */
1870 /* TRUE ... if equal (or dim == 1) */
1871 /* FALSE ... if unequal */
1872 /* */
1873 /* WARNING: */
1874 /* There is no checking of arguments in this function! */
1875 /*----------------------------------------------------------------------*/
1876 {
1877 /* There are (should be) only two possibilities:
1878 either all entries in the array point to the same distribution object;
1879 (set by unur_distr_cvec_set_marginals() call)
1880 or each entry has its own copy of some distribution object.
1881 (set by unur_distr_cvec_set_marginal_array() call)
1882 For dimension 1, TRUE is returned.
1883 */
1884 return (dim <= 1 || marginals[0] == marginals[1]) ? TRUE : FALSE;
1885 } /* end of _unur_distr_cvec_marginals_are_equal() */
1886
1887 /*---------------------------------------------------------------------------*/
1888
1889 int
_unur_distr_cvec_duplicate_firstmarginal(struct unur_distr * distr)1890 _unur_distr_cvec_duplicate_firstmarginal( struct unur_distr *distr )
1891 /*----------------------------------------------------------------------*/
1892 /* Duplicate first marginal distribution in array of marginal */
1893 /* distributions into all other slots of this array */
1894 /* This is only executed when all entries in this array point to the */
1895 /* same distribution object, i.e. when all marginal distributions */
1896 /* are equal. */
1897 /* */
1898 /* parameters: */
1899 /* distr ... pointer to distribution object */
1900 /* */
1901 /* return: */
1902 /* UNUR_SUCCESS ... on success */
1903 /* error code ... on error */
1904 /*----------------------------------------------------------------------*/
1905 {
1906 struct unur_distr *marginal;
1907 int i;
1908
1909 /* check arguments */
1910 CHECK_NULL( distr, UNUR_ERR_NULL );
1911 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1912
1913 marginal = DISTR.marginals[0];
1914
1915 /* marginal distributions known ? */
1916 if ( !(distr->set & UNUR_DISTR_SET_MARGINAL) || marginal==NULL ) {
1917 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"marginals");
1918 return UNUR_ERR_DISTR_DATA;
1919 }
1920
1921 /* marginal distribution are equal ? */
1922 if (!_unur_distr_cvec_marginals_are_equal(DISTR.marginals,distr->dim)) {
1923 /* nothing to do */
1924 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"marginals not equal");
1925 return UNUR_ERR_DISTR_DATA;
1926 }
1927
1928 /* make copy of marginal distribution objects */
1929 for (i=1; i<distr->dim; i++)
1930 DISTR.marginals[i] = _unur_distr_clone( marginal );
1931
1932 return UNUR_SUCCESS;
1933 } /* end of _unur_distr_cvec_duplicate_firstmarginal() */
1934
1935 /*---------------------------------------------------------------------------*/
1936
1937 int
unur_distr_cvec_set_pdfparams(struct unur_distr * distr,const double * params,int n_params)1938 unur_distr_cvec_set_pdfparams( struct unur_distr *distr, const double *params, int n_params )
1939 /*----------------------------------------------------------------------*/
1940 /* set array of parameters for distribution */
1941 /* */
1942 /* parameters: */
1943 /* distr ... pointer to distribution object */
1944 /* params ... list of arguments */
1945 /* n_params ... number of arguments */
1946 /* */
1947 /* return: */
1948 /* UNUR_SUCCESS ... on success */
1949 /* error code ... on error */
1950 /*----------------------------------------------------------------------*/
1951 {
1952 /* check arguments */
1953 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1954 _unur_check_NULL( NULL, params, UNUR_ERR_NULL );
1955 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
1956 if (n_params>0) _unur_check_NULL(distr->name,params,UNUR_ERR_NULL);
1957
1958 /* check number of new parameter for the distribution */
1959 if (n_params < 0 || n_params > UNUR_DISTR_MAXPARAMS ) {
1960 _unur_error(distr->name,UNUR_ERR_DISTR_NPARAMS,"");
1961 return UNUR_ERR_DISTR_NPARAMS;
1962 }
1963
1964 /* changelog */
1965 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
1966 /* derived parameters like mode, area, etc. might be wrong now! */
1967
1968 DISTR.n_params = n_params;
1969 if (n_params) memcpy( DISTR.params, params, n_params*sizeof(double) );
1970
1971 /* o.k. */
1972 return UNUR_SUCCESS;
1973 } /* end of unur_distr_cvec_set_pdfparams() */
1974
1975 /*---------------------------------------------------------------------------*/
1976
1977 int
unur_distr_cvec_get_pdfparams(const struct unur_distr * distr,const double ** params)1978 unur_distr_cvec_get_pdfparams( const struct unur_distr *distr, const double **params )
1979 /*----------------------------------------------------------------------*/
1980 /* get number of pdf parameters and sets pointer to array params[] of */
1981 /* parameters */
1982 /* */
1983 /* parameters: */
1984 /* distr ... pointer to distribution object */
1985 /* params ... pointer to list of arguments */
1986 /* */
1987 /* return: */
1988 /* number of pdf parameters */
1989 /* */
1990 /* error: */
1991 /* return 0 */
1992 /*----------------------------------------------------------------------*/
1993 {
1994 /* check arguments */
1995 _unur_check_NULL( NULL, distr, 0 );
1996 _unur_check_distr_object( distr, CVEC, 0 );
1997
1998 *params = (DISTR.n_params) ? DISTR.params : NULL;
1999 return DISTR.n_params;
2000
2001 } /* end of unur_distr_cvec_get_pdfparams() */
2002
2003
2004 /*---------------------------------------------------------------------------*/
2005
2006 int
unur_distr_cvec_set_pdfparams_vec(struct unur_distr * distr,int par,const double * param_vec,int n_param_vec)2007 unur_distr_cvec_set_pdfparams_vec( struct unur_distr *distr, int par, const double *param_vec, int n_param_vec )
2008 /*----------------------------------------------------------------------*/
2009 /* set vector array parameters for distribution */
2010 /* */
2011 /* parameters: */
2012 /* distr ... pointer to distribution object */
2013 /* par ... which parameter is set */
2014 /* param_vec ... parameter array with number `par' */
2015 /* n_param_vec ... length of parameter array */
2016 /* */
2017 /* return: */
2018 /* UNUR_SUCCESS ... on success */
2019 /* error code ... on error */
2020 /*----------------------------------------------------------------------*/
2021 {
2022 /* check arguments */
2023 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
2024 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
2025
2026 /* check new parameter for distribution */
2027 if (par < 0 || par >= UNUR_DISTR_MAXPARAMS ) {
2028 _unur_error(distr->name,UNUR_ERR_DISTR_NPARAMS,"");
2029 return UNUR_ERR_DISTR_NPARAMS;
2030 }
2031
2032
2033 if (param_vec != NULL) {
2034 /* allocate memory */
2035 DISTR.param_vecs[par] = _unur_xrealloc( DISTR.param_vecs[par], n_param_vec * sizeof(double) );
2036 /* copy parameters */
2037 memcpy( DISTR.param_vecs[par], param_vec, n_param_vec*sizeof(double) );
2038 /* set length of array */
2039 DISTR.n_param_vec[par] = n_param_vec;
2040 }
2041 else {
2042 if (DISTR.param_vecs[par]) free(DISTR.param_vecs[par]);
2043 DISTR.n_param_vec[par] = 0;
2044 }
2045
2046 /* changelog */
2047 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
2048 /* derived parameters like mode, area, etc. might be wrong now! */
2049
2050 /* o.k. */
2051 return UNUR_SUCCESS;
2052 } /* end of unur_distr_cvec_set_pdfparams_vec() */
2053
2054 /*---------------------------------------------------------------------------*/
2055
2056 int
unur_distr_cvec_get_pdfparams_vec(const struct unur_distr * distr,int par,const double ** param_vecs)2057 unur_distr_cvec_get_pdfparams_vec( const struct unur_distr *distr, int par, const double **param_vecs )
2058 /*----------------------------------------------------------------------*/
2059 /* get number of PDF parameters and sets pointer to array params[] of */
2060 /* parameters */
2061 /* */
2062 /* parameters: */
2063 /* distr ... pointer to distribution object */
2064 /* par ... which parameter is read */
2065 /* params ... pointer to parameter array with number `par' */
2066 /* */
2067 /* return: */
2068 /* length of parameter array with number `par' */
2069 /* */
2070 /* error: */
2071 /* return 0 */
2072 /*----------------------------------------------------------------------*/
2073 {
2074 /* check arguments */
2075 _unur_check_NULL( NULL, distr, 0 );
2076 _unur_check_distr_object( distr, CVEC, 0 );
2077
2078 /* check new parameter for distribution */
2079 if (par < 0 || par >= UNUR_DISTR_MAXPARAMS ) {
2080 _unur_error(distr->name,UNUR_ERR_DISTR_NPARAMS,"");
2081 *param_vecs = NULL;
2082 return 0;
2083 }
2084
2085 *param_vecs = DISTR.param_vecs[par];
2086
2087 return (*param_vecs) ? DISTR.n_param_vec[par] : 0;
2088 } /* end of unur_distr_cvec_get_pdfparams_vec() */
2089
2090 /*---------------------------------------------------------------------------*/
2091
2092 int
unur_distr_cvec_set_mode(struct unur_distr * distr,const double * mode)2093 unur_distr_cvec_set_mode( struct unur_distr *distr, const double *mode )
2094 /*----------------------------------------------------------------------*/
2095 /* set mode of distribution */
2096 /* */
2097 /* parameters: */
2098 /* distr ... pointer to distribution object */
2099 /* mode ... mode of PDF */
2100 /* */
2101 /* return: */
2102 /* UNUR_SUCCESS ... on success */
2103 /* error code ... on error */
2104 /*----------------------------------------------------------------------*/
2105 {
2106 int i;
2107
2108 /* check arguments */
2109 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
2110 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
2111
2112 /* we have to allocate memory first */
2113 if (DISTR.mode == NULL)
2114 DISTR.mode = _unur_xmalloc( distr->dim * sizeof(double) );
2115
2116 if (mode)
2117 /* mode vector given --> copy */
2118 memcpy( DISTR.mode, mode, distr->dim * sizeof(double) );
2119
2120 else /* mode == NULL --> use zero vector instead */
2121 for (i=0; i<distr->dim; i++)
2122 DISTR.mode[i] = 0.;
2123
2124 /* changelog */
2125 distr->set |= UNUR_DISTR_SET_MODE;
2126
2127 /* o.k. */
2128 return UNUR_SUCCESS;
2129 } /* end of unur_distr_cvec_set_mode() */
2130
2131 /*---------------------------------------------------------------------------*/
2132
2133 int
unur_distr_cvec_upd_mode(struct unur_distr * distr)2134 unur_distr_cvec_upd_mode( struct unur_distr *distr )
2135 /*----------------------------------------------------------------------*/
2136 /* (re-) compute mode of distribution (if possible) */
2137 /* */
2138 /* parameters: */
2139 /* distr ... pointer to distribution object */
2140 /* */
2141 /* return: */
2142 /* UNUR_SUCCESS ... on success */
2143 /* error code ... on error */
2144 /*----------------------------------------------------------------------*/
2145 {
2146 /* check arguments */
2147 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
2148 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
2149
2150 if (DISTR.upd_mode == NULL) {
2151 /* no function to compute mode available */
2152 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
2153 return UNUR_ERR_DISTR_DATA;
2154 }
2155
2156 /* compute mode */
2157 if ((DISTR.upd_mode)(distr)==UNUR_SUCCESS) {
2158 /* changelog */
2159 distr->set |= UNUR_DISTR_SET_MODE;
2160 return UNUR_SUCCESS;
2161 }
2162 else {
2163 /* computing of mode failed */
2164 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
2165 return UNUR_ERR_DISTR_DATA;
2166 }
2167
2168 } /* end of unur_distr_cvec_upd_mode() */
2169
2170 /*---------------------------------------------------------------------------*/
2171
2172 const double *
unur_distr_cvec_get_mode(struct unur_distr * distr)2173 unur_distr_cvec_get_mode( struct unur_distr *distr )
2174 /*----------------------------------------------------------------------*/
2175 /* get mode of distribution */
2176 /* */
2177 /* parameters: */
2178 /* distr ... pointer to distribution object */
2179 /* */
2180 /* return: */
2181 /* pointer to mode of distribution */
2182 /*----------------------------------------------------------------------*/
2183 {
2184 /* check arguments */
2185 _unur_check_NULL( NULL, distr, NULL );
2186 _unur_check_distr_object( distr, CVEC, NULL );
2187
2188 /* mode known ? */
2189 if ( !(distr->set & UNUR_DISTR_SET_MODE) ) {
2190 /* try to compute mode */
2191 if (DISTR.upd_mode == NULL) {
2192 /* no function to compute mode available */
2193 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"mode");
2194 return NULL;
2195 }
2196 else {
2197 /* compute mode */
2198 if (unur_distr_cvec_upd_mode(distr)!=UNUR_SUCCESS) {
2199 /* finding mode not successfully */
2200 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"mode");
2201 return NULL;
2202 }
2203 }
2204 }
2205
2206 return DISTR.mode;
2207
2208 } /* end of unur_distr_cvec_get_mode() */
2209
2210 /*---------------------------------------------------------------------------*/
2211
2212 int
unur_distr_cvec_set_center(struct unur_distr * distr,const double * center)2213 unur_distr_cvec_set_center( struct unur_distr *distr, const double *center )
2214 /*----------------------------------------------------------------------*/
2215 /* set center of distribution */
2216 /* */
2217 /* parameters: */
2218 /* distr ... pointer to distribution object */
2219 /* center ... center of PDF */
2220 /* */
2221 /* return: */
2222 /* UNUR_SUCCESS ... on success */
2223 /* error code ... on error */
2224 /*----------------------------------------------------------------------*/
2225 {
2226 int i;
2227
2228 /* check arguments */
2229 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
2230 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
2231
2232 /* we have to allocate memory first */
2233 if (DISTR.center == NULL)
2234 DISTR.center = _unur_xmalloc( distr->dim * sizeof(double) );
2235
2236 if (center)
2237 /* center vector given --> copy */
2238 memcpy( DISTR.center, center, distr->dim * sizeof(double) );
2239
2240 else /* center == NULL --> use zero vector instead */
2241 for (i=0; i<distr->dim; i++)
2242 DISTR.center[i] = 0.;
2243
2244 /* changelog */
2245 distr->set |= UNUR_DISTR_SET_CENTER;
2246
2247 /* o.k. */
2248 return UNUR_SUCCESS;
2249 } /* end of unur_distr_cvec_set_center() */
2250
2251 /*---------------------------------------------------------------------------*/
2252
2253 const double *
unur_distr_cvec_get_center(struct unur_distr * distr)2254 unur_distr_cvec_get_center( struct unur_distr *distr )
2255 /*----------------------------------------------------------------------*/
2256 /* get center of distribution */
2257 /* */
2258 /* parameters: */
2259 /* distr ... pointer to distribution object */
2260 /* */
2261 /* return: */
2262 /* pointer to center of distribution */
2263 /*----------------------------------------------------------------------*/
2264 {
2265 int i;
2266
2267 /* check arguments */
2268 _unur_check_NULL( NULL, distr, NULL );
2269 _unur_check_distr_object( distr, CVEC, NULL );
2270
2271 /* center given */
2272 if ( distr->set & UNUR_DISTR_SET_CENTER )
2273 return DISTR.center;
2274
2275 /* else try mode */
2276 if ( distr->set & UNUR_DISTR_SET_MODE )
2277 return DISTR.mode;
2278
2279 /* else try mean */
2280 if ( distr->set & UNUR_DISTR_SET_MEAN )
2281 return DISTR.mean;
2282
2283 /* otherwise use (0,...,0) */
2284 if ( DISTR.center == NULL )
2285 DISTR.center = _unur_xmalloc( distr->dim * sizeof(double) );
2286 for (i=0; i<distr->dim; i++)
2287 DISTR.center[i] = 0.;
2288
2289 return DISTR.center;
2290
2291 } /* end of unur_distr_cvec_get_center() */
2292
2293 /*---------------------------------------------------------------------------*/
2294
2295 int
unur_distr_cvec_set_pdfvol(struct unur_distr * distr,double volume)2296 unur_distr_cvec_set_pdfvol( struct unur_distr *distr, double volume )
2297 /*----------------------------------------------------------------------*/
2298 /* set volume below PDF */
2299 /* */
2300 /* parameters: */
2301 /* distr ... pointer to distribution object */
2302 /* volume ... volume below PDF */
2303 /* */
2304 /* return: */
2305 /* UNUR_SUCCESS ... on success */
2306 /* error code ... on error */
2307 /*----------------------------------------------------------------------*/
2308 {
2309 /* check arguments */
2310 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
2311 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
2312
2313 /* check new parameter for distribution */
2314 if (volume <= 0.) {
2315 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"PDF volume <= 0");
2316 return UNUR_ERR_DISTR_SET;
2317 }
2318
2319 DISTR.volume = volume;
2320
2321 /* changelog */
2322 distr->set |= UNUR_DISTR_SET_PDFVOLUME;
2323
2324 /* o.k. */
2325 return UNUR_SUCCESS;
2326
2327 } /* end of unur_distr_cvec_set_pdfvol() */
2328
2329 /*---------------------------------------------------------------------------*/
2330
2331 int
unur_distr_cvec_upd_pdfvol(struct unur_distr * distr)2332 unur_distr_cvec_upd_pdfvol( struct unur_distr *distr )
2333 /*----------------------------------------------------------------------*/
2334 /* (re-) compute volume below p.d.f. of distribution (if possible) */
2335 /* */
2336 /* parameters: */
2337 /* distr ... pointer to distribution object */
2338 /* */
2339 /* return: */
2340 /* UNUR_SUCCESS ... on success */
2341 /* error code ... on error */
2342 /*----------------------------------------------------------------------*/
2343 {
2344 /* check arguments */
2345 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
2346 _unur_check_distr_object( distr, CVEC, UNUR_ERR_DISTR_INVALID );
2347
2348 if (DISTR.upd_volume == NULL) {
2349 /* no function to compute mode available */
2350 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
2351 return UNUR_ERR_DISTR_DATA;
2352 }
2353
2354 /* compute area */
2355 if (((DISTR.upd_volume)(distr)!=UNUR_SUCCESS) || DISTR.volume <= 0.) {
2356 /* computing of area failed */
2357 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"upd volume <= 0");
2358 DISTR.volume = 1.; /* avoid possible floating point exceptions */
2359 distr->set &= ~UNUR_DISTR_SET_PDFVOLUME;
2360 return UNUR_ERR_DISTR_SET;
2361 }
2362
2363 /* changelog */
2364 distr->set |= UNUR_DISTR_SET_PDFVOLUME;
2365
2366 return UNUR_SUCCESS;
2367 } /* end of unur_distr_cvec_upd_pdfvol() */
2368
2369 /*---------------------------------------------------------------------------*/
2370
2371 double
unur_distr_cvec_get_pdfvol(struct unur_distr * distr)2372 unur_distr_cvec_get_pdfvol( struct unur_distr *distr )
2373 /*----------------------------------------------------------------------*/
2374 /* get volume below PDF of distribution */
2375 /* */
2376 /* parameters: */
2377 /* distr ... pointer to distribution object */
2378 /* */
2379 /* return: */
2380 /* volume below PDF of distribution */
2381 /*----------------------------------------------------------------------*/
2382 {
2383 /* check arguments */
2384 _unur_check_NULL( NULL, distr, INFINITY );
2385 _unur_check_distr_object( distr, CVEC, INFINITY );
2386
2387 /* volume known ? */
2388 if ( !(distr->set & UNUR_DISTR_SET_PDFVOLUME) ) {
2389 /* try to compute volume */
2390 if (DISTR.upd_volume == NULL) {
2391 /* no function to compute area available */
2392 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"volume");
2393 return INFINITY;
2394 }
2395 else {
2396 /* compute area */
2397 unur_distr_cvec_upd_pdfvol( distr );
2398 }
2399 }
2400
2401 return DISTR.volume;
2402
2403 } /* end of unur_distr_cvec_get_pdfvol() */
2404
2405 /*****************************************************************************/
2406 /* call PDFs and their derivatives */
2407 /* (internal functions. no checking for NULL pointer) */
2408
2409 double
_unur_cvec_PDF(const double * x,struct unur_distr * distr)2410 _unur_cvec_PDF(const double *x, struct unur_distr *distr)
2411 {
2412 if ( (distr->set & UNUR_DISTR_SET_DOMAINBOUNDED) &&
2413 !_unur_distr_cvec_is_indomain(x, distr) )
2414 return 0.;
2415
2416 return (*(distr->data.cvec.pdf)) (x,distr);
2417 }
2418
2419 int
_unur_cvec_dPDF(double * result,const double * x,struct unur_distr * distr)2420 _unur_cvec_dPDF(double *result, const double *x, struct unur_distr *distr)
2421 {
2422 int d;
2423 if ( (distr->set & UNUR_DISTR_SET_DOMAINBOUNDED) &&
2424 !_unur_distr_cvec_is_indomain(x, distr) ) {
2425 for (d=0; d < distr->dim; d++) result[d] = 0.;
2426 return UNUR_SUCCESS;
2427 }
2428
2429 return (*(distr->data.cvec.dpdf)) (result,x,distr);
2430 }
2431
2432 double
_unur_cvec_pdPDF(const double * x,int coord,struct unur_distr * distr)2433 _unur_cvec_pdPDF(const double *x, int coord, struct unur_distr *distr)
2434 {
2435 if ( (distr->set & UNUR_DISTR_SET_DOMAINBOUNDED) &&
2436 !_unur_distr_cvec_is_indomain(x, distr) )
2437 return 0.;
2438
2439 return (*(distr->data.cvec.pdpdf)) (x,coord,distr);
2440 }
2441
2442 double
_unur_cvec_logPDF(const double * x,struct unur_distr * distr)2443 _unur_cvec_logPDF(const double *x, struct unur_distr *distr)
2444 {
2445 if ( (distr->set & UNUR_DISTR_SET_DOMAINBOUNDED) &&
2446 !_unur_distr_cvec_is_indomain(x, distr) )
2447 return -INFINITY;
2448
2449 return (*(distr->data.cvec.logpdf)) (x,distr);
2450 }
2451
2452 int
_unur_cvec_dlogPDF(double * result,const double * x,struct unur_distr * distr)2453 _unur_cvec_dlogPDF(double *result, const double *x, struct unur_distr *distr)
2454 {
2455 int d;
2456 if ( (distr->set & UNUR_DISTR_SET_DOMAINBOUNDED) &&
2457 !_unur_distr_cvec_is_indomain(x, distr) ) {
2458 for (d=0; d < distr->dim; d++) result[d] = 0.;
2459 return UNUR_SUCCESS;
2460 }
2461
2462 return (*(distr->data.cvec.dlogpdf)) (result,x,distr);
2463 }
2464
2465 double
_unur_cvec_pdlogPDF(const double * x,int coord,struct unur_distr * distr)2466 _unur_cvec_pdlogPDF(const double *x, int coord, struct unur_distr *distr)
2467 {
2468 if ( (distr->set & UNUR_DISTR_SET_DOMAINBOUNDED) &&
2469 !_unur_distr_cvec_is_indomain(x, distr) )
2470 return 0.;
2471
2472 return (*(distr->data.cvec.pdlogpdf)) (x,coord,distr);
2473 }
2474
2475 /*****************************************************************************/
2476
2477 /*---------------------------------------------------------------------------*/
2478 #ifdef UNUR_ENABLE_LOGGING
2479 /*---------------------------------------------------------------------------*/
2480
2481 void
_unur_distr_cvec_debug(const struct unur_distr * distr,const char * genid)2482 _unur_distr_cvec_debug( const struct unur_distr *distr, const char *genid )
2483 /*----------------------------------------------------------------------*/
2484 /* write info about distribution into LOG file */
2485 /* */
2486 /* parameters: */
2487 /* distr ... pointer to distribution object */
2488 /* genid ... pointer to generator id */
2489 /*----------------------------------------------------------------------*/
2490 {
2491 FILE *LOG;
2492 double *mat;
2493
2494 /* check arguments */
2495 CHECK_NULL(distr,RETURN_VOID);
2496 COOKIE_CHECK(distr,CK_DISTR_CVEC,RETURN_VOID);
2497
2498 LOG = unur_get_stream();
2499
2500 fprintf(LOG,"%s: distribution:\n",genid);
2501 fprintf(LOG,"%s:\ttype = continuous multivariate distribution\n",genid);
2502 fprintf(LOG,"%s:\tname = %s\n",genid,distr->name);
2503
2504 fprintf(LOG,"%s:\tdimension = %d\n",genid,distr->dim);
2505
2506 fprintf(LOG,"%s:\tfunctions: ",genid);
2507 if (DISTR.pdf) fprintf(LOG,"PDF ");
2508 if (DISTR.logpdf) fprintf(LOG,"logPDF ");
2509 if (DISTR.dpdf) fprintf(LOG,"dPDF ");
2510 if (DISTR.dlogpdf) fprintf(LOG,"dlogPDF ");
2511 if (DISTR.pdpdf) fprintf(LOG,"pdPDF ");
2512 if (DISTR.pdlogpdf) fprintf(LOG,"pdlogPDF ");
2513 fprintf(LOG,"\n%s:\n",genid);
2514
2515 /* domain */
2516 fprintf(LOG,"%s:\tdomain = ",genid);
2517 if (!(distr->set & UNUR_DISTR_SET_DOMAINBOUNDED)) {
2518 fprintf(LOG,"unbounded\n");
2519 }
2520 else {
2521 if (DISTR.domainrect) {
2522 double *domain = DISTR.domainrect;
2523 int i;
2524 fprintf(LOG,"rectangular\n");
2525 for (i=0; i<distr->dim; i++)
2526 fprintf(LOG,"%s:\t %c ( %g, %g)\n",genid, i?'x':' ',
2527 domain[2*i], domain[2*i+1]);
2528 }
2529 }
2530 fprintf(LOG,"%s:\n",genid);
2531
2532 /* mode */
2533 mat = ((distr->set & UNUR_DISTR_SET_MODE) && DISTR.mode) ? DISTR.mode : NULL;
2534 _unur_matrix_print_vector( distr->dim, mat, "\tmode =", LOG, genid, "\t ");
2535
2536 /* mean vector */
2537 mat = ((distr->set & UNUR_DISTR_SET_MEAN) && DISTR.mean) ? DISTR.mean : NULL;
2538 _unur_matrix_print_vector( distr->dim, mat, "\tmean vector =", LOG, genid, "\t ");
2539
2540 /* center vector */
2541 if ((distr->set & UNUR_DISTR_SET_CENTER) && DISTR.center)
2542 _unur_matrix_print_vector( distr->dim, DISTR.center, "\tcenter vector =", LOG, genid, "\t ");
2543 else {
2544 fprintf(LOG,"%s:\tcenter = mode [not given explicitly]\n",genid);
2545 fprintf(LOG,"%s:\n",genid);
2546 }
2547
2548 /* covariance matrix */
2549 mat = ((distr->set & UNUR_DISTR_SET_COVAR) && DISTR.covar) ? DISTR.covar : NULL;
2550 _unur_matrix_print_matrix( distr->dim, mat, "\tcovariance matrix =", LOG, genid, "\t ");
2551
2552 /* inverse covariance matrix */
2553 /* mat = ((distr->set & UNUR_DISTR_SET_COVAR_INV) && DISTR.covar_inv) ? DISTR.covar_inv : NULL; */
2554 /* _unur_matrix_print_matrix( distr->dim, mat, "\tinverse covariance matrix =", LOG, genid, "\t "); */
2555
2556 /* cholesky factor of covariance matrix */
2557 mat = ((distr->set & UNUR_DISTR_SET_CHOLESKY) && DISTR.cholesky) ? DISTR.cholesky : NULL;
2558 _unur_matrix_print_matrix( distr->dim, mat, "\tcholesky factor of covariance matrix =", LOG, genid, "\t ");
2559
2560 /* rank correlation matrix */
2561 mat = ((distr->set & UNUR_DISTR_SET_RANKCORR) && DISTR.rankcorr) ? DISTR.rankcorr : NULL;
2562 _unur_matrix_print_matrix( distr->dim, mat, "\trank correlation matrix =", LOG, genid, "\t ");
2563
2564 /* cholesky factor of rank correlation matrix */
2565 mat = ((distr->set & UNUR_DISTR_SET_RK_CHOLESKY) && DISTR.rk_cholesky) ? DISTR.rk_cholesky : NULL;
2566 _unur_matrix_print_matrix( distr->dim, mat, "\tcholesky factor of rank correlation matrix =", LOG, genid, "\t ");
2567
2568 /* marginal distributions */
2569 fprintf(LOG,"%s:\tmarginal distributions:\n",genid);
2570 if (distr->set & UNUR_DISTR_SET_MARGINAL) {
2571 if (_unur_distr_cvec_marginals_are_equal(DISTR.marginals, distr->dim)) {
2572 fprintf(LOG,"%s: all mariginals [1-%d]:\n",genid,distr->dim);
2573 _unur_distr_cont_debug( DISTR.marginals[0], genid );
2574 }
2575 else {
2576 int i;
2577 for (i=0; i<distr->dim; i++) {
2578 fprintf(LOG,"%s: mariginal [%d]:\n",genid,i+1);
2579 _unur_distr_cont_debug( DISTR.marginals[i], genid );
2580 }
2581 }
2582 }
2583 else {
2584 fprintf(LOG,"%s:\t [unknown]\n",genid);
2585 }
2586 fprintf(LOG,"%s:\n",genid);
2587
2588 #ifdef USE_DEPRECATED_CODE
2589 /* standardized marginal distributions */
2590 fprintf(LOG,"%s:\tstandardized marginal distributions: [see also marginal distributions]\n",genid);
2591 if (distr->set & UNUR_DISTR_SET_STDMARGINAL) {
2592 if (_unur_distr_cvec_marginals_are_equal(DISTR.stdmarginals, distr->dim)) {
2593 fprintf(LOG,"%s: all standardized mariginals [1-%d]:\n",genid,distr->dim);
2594 _unur_distr_cont_debug( DISTR.stdmarginals[0], genid );
2595 }
2596 else {
2597 int i;
2598 for (i=0; i<distr->dim; i++) {
2599 fprintf(LOG,"%s: mariginal [%d]:\n",genid,i+1);
2600 _unur_distr_cont_debug( DISTR.stdmarginals[i], genid );
2601 }
2602 }
2603 }
2604 else {
2605 fprintf(LOG,"%s:\t [unknown]\n",genid);
2606 }
2607 fprintf(LOG,"%s:\n",genid);
2608 #endif
2609
2610 } /* end of _unur_distr_cvec_debug() */
2611
2612 /*---------------------------------------------------------------------------*/
2613 #endif /* end UNUR_ENABLE_LOGGING */
2614 /*---------------------------------------------------------------------------*/
2615
2616 /*---------------------------------------------------------------------------*/
2617 #undef DISTR
2618 /*---------------------------------------------------------------------------*/
2619
2620