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