1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: discr.c *
8 * *
9 * manipulate univariate discrete 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 "distr_source.h"
43 #include "distr.h"
44 #include "discr.h"
45
46 /*---------------------------------------------------------------------------*/
47
48 #define DISTR distr->data.discr
49
50 /*---------------------------------------------------------------------------*/
51 /* Constants */
52
53 /* maximum size of domain for which the pmfsum is computed automatically */
54 #define MAX_PMF_DOMAIN_FOR_UPD_PMFSUM (1000)
55
56 /*---------------------------------------------------------------------------*/
57
58 static double _unur_distr_discr_eval_pmf_tree( int k, const struct unur_distr *distr );
59 /*---------------------------------------------------------------------------*/
60 /* evaluate function tree for PMF. */
61 /*---------------------------------------------------------------------------*/
62
63 static double _unur_distr_discr_eval_cdf_tree( int k, const struct unur_distr *distr );
64 /*---------------------------------------------------------------------------*/
65 /* evaluate function tree for CDF. */
66 /*---------------------------------------------------------------------------*/
67
68 static void _unur_distr_discr_free( struct unur_distr *distr );
69 /*---------------------------------------------------------------------------*/
70 /* destroy distribution object. */
71 /*---------------------------------------------------------------------------*/
72
73 /*---------------------------------------------------------------------------*/
74
75 static int _unur_distr_discr_find_mode( struct unur_distr *distr );
76 /*---------------------------------------------------------------------------*/
77 /* find mode of unimodal probability vector numerically by bisection. */
78 /*---------------------------------------------------------------------------*/
79
80 /*****************************************************************************/
81 /** **/
82 /** univariate discrete distributions **/
83 /** **/
84 /*****************************************************************************/
85
86 /*---------------------------------------------------------------------------*/
87
88 struct unur_distr *
unur_distr_discr_new(void)89 unur_distr_discr_new( void )
90 /*----------------------------------------------------------------------*/
91 /* create a new (empty) distribution object */
92 /* type: univariate discete */
93 /* */
94 /* parameters: */
95 /* none */
96 /* */
97 /* return: */
98 /* pointer to distribution object */
99 /* */
100 /* error: */
101 /* return NULL */
102 /*----------------------------------------------------------------------*/
103 {
104 register struct unur_distr *distr;
105 register int i;
106
107 /* get empty distribution object */
108 distr = _unur_distr_generic_new();
109 if (!distr) return NULL;
110
111 /* set magic cookie */
112 COOKIE_SET(distr,CK_DISTR_DISCR);
113
114 /* set type of distribution */
115 distr->type = UNUR_DISTR_DISCR;
116
117 /* set id to generic distribution */
118 distr->id = UNUR_DISTR_GENERIC;
119
120 /* dimension of random vector */
121 distr->dim = 1; /* univariant */
122
123 /* destructor */
124 distr->destroy = _unur_distr_discr_free;
125
126 /* clone */
127 distr->clone = _unur_distr_discr_clone;
128
129 /* set defaults */
130
131 /* finite probability vector */
132 DISTR.pv = NULL; /* probability vector (PV) */
133 DISTR.n_pv = 0; /* length of PV */
134
135 /* probability mass function */
136 DISTR.pmf = NULL; /* pointer to PMF */
137 DISTR.cdf = NULL; /* pointer to CDF */
138 DISTR.invcdf = NULL; /* pointer to inverse CDF */
139
140 DISTR.init = NULL; /* pointer to special init routine */
141
142 DISTR.set_params= NULL; /* funct for setting parameters and domain*/
143
144 DISTR.n_params = 0; /* number of parameters of the pmf */
145 /* initialize parameters of the PMF */
146 for (i=0; i<UNUR_DISTR_MAXPARAMS; i++)
147 DISTR.params[i] = 0.;
148
149 DISTR.norm_constant = 1.; /* (log of) normalization constant for PMF
150 (initialized to avoid accidently floating
151 point exception */
152
153 DISTR.trunc[0] = DISTR.domain[0] = 0; /* left boundary of domain */
154 DISTR.trunc[1] = DISTR.domain[1] = INT_MAX; /* right boundary of domain */
155
156 DISTR.mode = 0; /* location of mode */
157 DISTR.upd_mode = _unur_distr_discr_find_mode; /* funct for computing mode */
158
159 DISTR.sum = 1.; /* sum over PMF */
160 DISTR.upd_sum = NULL; /* funct for computing sum */
161
162 DISTR.pmftree = NULL; /* pointer to function tree for PMF */
163 DISTR.cdftree = NULL; /* pointer to function tree for CDF */
164
165 /* return pointer to object */
166 return distr;
167
168 } /* end of unur_distr_discr_new() */
169
170 /*---------------------------------------------------------------------------*/
171
172 struct unur_distr *
_unur_distr_discr_clone(const struct unur_distr * distr)173 _unur_distr_discr_clone( const struct unur_distr *distr )
174 /*----------------------------------------------------------------------*/
175 /* copy (clone) distribution object */
176 /* */
177 /* parameters: */
178 /* distr ... pointer to source distribution object */
179 /* */
180 /* return: */
181 /* pointer to clone of distribution object */
182 /*----------------------------------------------------------------------*/
183 {
184 #define CLONE clone->data.discr
185
186 struct unur_distr *clone;
187
188 /* check arguments */
189 _unur_check_NULL( NULL, distr, NULL );
190 _unur_check_distr_object( distr, DISCR, NULL );
191
192 /* allocate memory */
193 clone = _unur_xmalloc( sizeof(struct unur_distr) );
194
195 /* copy distribution object into clone */
196 memcpy( clone, distr, sizeof( struct unur_distr ) );
197
198 /* copy function trees into generator object (when there is one) */
199 CLONE.pmftree = (DISTR.pmftree) ? _unur_fstr_dup_tree(DISTR.pmftree) : NULL;
200 CLONE.cdftree = (DISTR.cdftree) ? _unur_fstr_dup_tree(DISTR.cdftree) : NULL;
201
202 /* copy probability vector into generator object (when there is one) */
203 if (DISTR.pv) {
204 CLONE.pv = _unur_xmalloc( DISTR.n_pv * sizeof(double) );
205 memcpy( CLONE.pv, DISTR.pv, DISTR.n_pv * sizeof(double) );
206 }
207
208 /* copy user name for distribution */
209 if (distr->name_str) {
210 size_t len = strlen(distr->name_str) + 1;
211 clone->name_str = _unur_xmalloc(len);
212 memcpy( clone->name_str, distr->name_str, len );
213 clone->name = clone->name_str;
214 }
215
216 return clone;
217
218 #undef CLONE
219 } /* end of _unur_distr_discr_clone() */
220
221 /*---------------------------------------------------------------------------*/
222
223 void
_unur_distr_discr_free(struct unur_distr * distr)224 _unur_distr_discr_free( struct unur_distr *distr )
225 /*----------------------------------------------------------------------*/
226 /* free distribution object */
227 /* */
228 /* parameters: */
229 /* distr ... pointer to distribution object */
230 /*----------------------------------------------------------------------*/
231 {
232 /* check arguments */
233 if( distr == NULL ) /* nothing to do */
234 return;
235 _unur_check_distr_object( distr, DISCR, RETURN_VOID );
236
237 if (DISTR.pmftree) _unur_fstr_free(DISTR.pmftree);
238 if (DISTR.cdftree) _unur_fstr_free(DISTR.cdftree);
239
240 if (DISTR.pv) free( DISTR.pv );
241
242 /* user name for distribution */
243 if (distr->name_str) free(distr->name_str);
244
245 COOKIE_CLEAR(distr);
246 free( distr );
247
248 } /* end of unur_distr_discr_free() */
249
250 /*---------------------------------------------------------------------------*/
251
252 int
unur_distr_discr_set_pv(struct unur_distr * distr,const double * pv,int n_pv)253 unur_distr_discr_set_pv( struct unur_distr *distr, const double *pv, int n_pv )
254 /*----------------------------------------------------------------------*/
255 /* set probability vector for distribution */
256 /* */
257 /* parameters: */
258 /* distr ... pointer to distribution object */
259 /* pv ... pointer to PV */
260 /* n_pv ... length of PV */
261 /* */
262 /* return: */
263 /* UNUR_SUCCESS ... on success */
264 /* error code ... on error */
265 /*----------------------------------------------------------------------*/
266 {
267 /* check arguments */
268 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
269 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
270
271 /* it is not possible to set a PV when a PMF is given. */
272 if (DISTR.pmf != NULL || DISTR.cdf != NULL) {
273 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"PMF/CDF given, cannot set PV");
274 return UNUR_ERR_DISTR_SET;
275 }
276
277 /* check new parameter for distribution */
278 if (n_pv < 0) {
279 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"length of PV");
280 return UNUR_ERR_DISTR_SET;
281 }
282
283 /* n_pv must not be too large */
284 if ( (DISTR.domain[0] > 0) && ((unsigned)DISTR.domain[0] + (unsigned)n_pv > INT_MAX) ) {
285 /* n_pv too large, causes overflow */
286 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"length of PV too large, overflow");
287 return UNUR_ERR_DISTR_SET;
288 }
289 DISTR.domain[1] = DISTR.domain[0] + n_pv - 1;
290
291 /* we do not check non-negativity of p.v.
292 (it is cheaper to do it when unur_init() is called */
293
294 /* allocate memory for probability vector */
295 DISTR.pv = _unur_xrealloc( DISTR.pv, n_pv * sizeof(double) );
296 if (!DISTR.pv) return UNUR_ERR_MALLOC;
297
298 /* copy probability vector */
299 memcpy( DISTR.pv, pv, n_pv * sizeof(double) );
300 DISTR.n_pv = n_pv;
301
302 /* o.k. */
303 return UNUR_SUCCESS;
304 } /* end of unur_distr_discr_set_pv() */
305
306 /*---------------------------------------------------------------------------*/
307
308 int
unur_distr_discr_make_pv(struct unur_distr * distr)309 unur_distr_discr_make_pv( struct unur_distr *distr )
310 /*----------------------------------------------------------------------*/
311 /* compute probability vector */
312 /* */
313 /* parameters: */
314 /* distr ... pointer to distribution object */
315 /* */
316 /* return: */
317 /* length of probability vector */
318 /* */
319 /* error: */
320 /* return 0 */
321 /*----------------------------------------------------------------------*/
322 {
323 double *pv; /* pointer to probability vector */
324 int n_pv; /* length of PV */
325 double cdf, cdf_old; /* cumulated sum of PV */
326 double thresh_cdf; /* threshold for truncating PV */
327 int valid; /* whether cumputed PV is valid */
328 int i;
329
330 /* check arguments */
331 _unur_check_NULL( NULL, distr, 0 );
332 _unur_check_distr_object( distr, DISCR, 0 );
333
334 /* PMF or CDF required */
335 if ( DISTR.pmf == NULL && DISTR.cdf == NULL) {
336 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"PMF or CDF");
337 return 0;
338 }
339
340 /* if there exists a PV, it has to be removed */
341 if (DISTR.pv != NULL) {
342 free(DISTR.pv); DISTR.n_pv = 0;
343 }
344
345 /* compute PV */
346
347 if ((unsigned)DISTR.domain[1] - (unsigned)DISTR.domain[0] < UNUR_MAX_AUTO_PV ) {
348
349 /* first case: bounded domain */
350 n_pv = DISTR.domain[1] - DISTR.domain[0] + 1;
351 pv = _unur_xmalloc( n_pv * sizeof(double) );
352 if (DISTR.pmf) {
353 for (i=0; i<n_pv; i++)
354 pv[i] = _unur_discr_PMF(DISTR.domain[0]+i,distr);
355 }
356 else if (DISTR.cdf) {
357 cdf_old = 0.;
358 for (i=0; i<n_pv; i++) {
359 cdf = _unur_discr_CDF(DISTR.domain[0]+i,distr);
360 pv[i] = cdf - cdf_old;
361 cdf_old = cdf;
362 }
363 }
364 valid = TRUE;
365 }
366
367 else {
368 /* second case: domain too big but sum over PMF given */
369 /* we chop off the trailing part of the distribution */
370 #define MALLOC_SIZE 1000 /* allocate 1000 doubles at once */
371
372 int n_alloc; /* number of doubles allocated */
373 int max_alloc; /* maximal number of allocated doubles */
374 int size_alloc; /* size of allocated blocks */
375
376 /* get maximal size of PV */
377 if ( (DISTR.domain[0] <= 0) || (INT_MAX - DISTR.domain[0] >= UNUR_MAX_AUTO_PV - 1) ) {
378 /* we can have a PV of length UNUR_MAX_AUTO_PV */
379 size_alloc = MALLOC_SIZE;
380 max_alloc = UNUR_MAX_AUTO_PV;
381 }
382 else { /* length of PV must be shorter than UNUR_MAX_AUTO_PV */
383 size_alloc = max_alloc = INT_MAX - DISTR.domain[0];
384 }
385
386 /* init counter */
387 n_pv = 0;
388 pv = NULL;
389 valid = FALSE; /* created PV is empty yet and not valid */
390 cdf = 0.; /* cumulated sum of PV */
391 cdf_old = 0.; /* cumulated sum of PV in last iteration */
392 /* threshold for truncating PV */
393 thresh_cdf = (distr->set & UNUR_DISTR_SET_PMFSUM) ? (1.-1.e-8)*DISTR.sum : INFINITY;
394
395 /* compute PV */
396 for (n_alloc = size_alloc; n_alloc <= max_alloc; n_alloc += size_alloc) {
397 pv = _unur_xrealloc( pv, n_alloc * sizeof(double) );
398
399 if (DISTR.pmf) {
400 for (i=0; i<size_alloc; i++) {
401 cdf += pv[n_pv] = _unur_discr_PMF(DISTR.domain[0]+n_pv,distr);
402 n_pv++;
403 if (cdf > thresh_cdf) { valid = TRUE; break; }
404 }
405 }
406 else if (DISTR.cdf) {
407 for (i=0; i<size_alloc; i++) {
408 cdf = _unur_discr_CDF(DISTR.domain[0]+n_pv,distr);
409 pv[n_pv] = cdf - cdf_old;
410 cdf_old = cdf;
411 n_pv++;
412 if (cdf > thresh_cdf) { valid = TRUE; break; }
413 }
414 }
415 if (cdf > thresh_cdf) break;
416 }
417
418 if (distr->set & UNUR_DISTR_SET_PMFSUM) {
419 /* make a warning if computed PV might not be valid */
420 if (valid != TRUE)
421 /* not successful */
422 _unur_warning(distr->name,UNUR_ERR_DISTR_GET,"PV truncated");
423 }
424 else { /* PMFSUM not known */
425 /* assume we have the important part of distribution */
426 valid = TRUE;
427 DISTR.sum = cdf;
428 distr->set |= UNUR_DISTR_SET_PMFSUM;
429 }
430
431 #undef MALLOC_SIZE
432 }
433
434 /* store vector */
435 DISTR.pv = pv;
436 DISTR.n_pv = n_pv;
437 DISTR.domain[1] = DISTR.domain[0] + n_pv - 1;
438
439 /* o.k. */
440 return (valid) ? n_pv : -n_pv;
441 } /* end of unur_distr_discr_make_pv() */
442
443 /*---------------------------------------------------------------------------*/
444
445 int
unur_distr_discr_get_pv(const struct unur_distr * distr,const double ** pv)446 unur_distr_discr_get_pv( const struct unur_distr *distr, const double **pv )
447 /*----------------------------------------------------------------------*/
448 /* get length of probability vector and set pointer to probability */
449 /* vector */
450 /* */
451 /* parameters: */
452 /* distr ... pointer to distribution object */
453 /* pv ... pointer to probability vector */
454 /* */
455 /* return: */
456 /* length of probability vector */
457 /* */
458 /* error: */
459 /* return 0 */
460 /*----------------------------------------------------------------------*/
461 {
462 /* check arguments */
463 _unur_check_NULL( NULL, distr, 0 );
464 _unur_check_distr_object( distr, DISCR, 0 );
465
466 *pv = (DISTR.pv) ? DISTR.pv : NULL;
467 return DISTR.n_pv;
468
469 } /* end of unur_distr_discr_get_pv() */
470
471 /*---------------------------------------------------------------------------*/
472
473 double
unur_distr_discr_eval_pv(int k,const struct unur_distr * distr)474 unur_distr_discr_eval_pv( int k, const struct unur_distr *distr )
475 /*----------------------------------------------------------------------*/
476 /* returns the value of the probability vector at k or, if there is no */
477 /* probability vector defined, evaluates the pmf */
478 /* */
479 /* parampeters: */
480 /* k ... argument for probability vector of pmf */
481 /* distr ... pointer to distribution object */
482 /* */
483 /* return: */
484 /* pv[k] or pmf(k) */
485 /*----------------------------------------------------------------------*/
486 {
487 /* check arguments */
488 _unur_check_NULL( NULL, distr, INFINITY );
489 _unur_check_distr_object( distr, DISCR, INFINITY );
490
491 if (DISTR.pv != NULL) {
492 /* use probability vector */
493 if (k < DISTR.domain[0] || k > DISTR.domain[1])
494 return 0.;
495 else
496 return (DISTR.pv[k-DISTR.domain[0]]);
497 }
498
499 if (DISTR.pmf != NULL) {
500 /* use PMF */
501 double px = _unur_discr_PMF(k,distr);
502 if (_unur_isnan(px)) {
503 _unur_warning(distr->name,UNUR_ERR_DISTR_DATA,"PMF returns NaN");
504 return 0.;
505 }
506 else
507 return px;
508 }
509
510 /* else: data missing */
511 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
512 return INFINITY;
513
514 } /* end of unur_distr_discr_eval_pv() */
515
516 /*---------------------------------------------------------------------------*/
517
518 int
unur_distr_discr_set_pmf(struct unur_distr * distr,UNUR_FUNCT_DISCR * pmf)519 unur_distr_discr_set_pmf( struct unur_distr *distr, UNUR_FUNCT_DISCR *pmf )
520 /*----------------------------------------------------------------------*/
521 /* set PMF of distribution */
522 /* */
523 /* parameters: */
524 /* distr ... pointer to distribution object */
525 /* pmf ... pointer to PMF */
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, pmf, UNUR_ERR_NULL );
535 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
536
537 /* it is not possible to set both a PMF and a PV */
538 if (DISTR.pv != NULL) {
539 _unur_warning(distr->name,UNUR_ERR_DISTR_SET,"delete exisiting PV");
540 free(DISTR.pv); DISTR.n_pv = 0;
541 }
542
543 /* we do not allow overwriting a PMF */
544 if (DISTR.pmf != NULL) {
545 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of PMF not allowed");
546 return UNUR_ERR_DISTR_SET;
547 }
548
549 /* changelog */
550 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
551 /* derived parameters like mode, sum, etc. might be wrong now! */
552
553 DISTR.pmf = pmf;
554 return UNUR_SUCCESS;
555
556 } /* end of unur_distr_discr_set_pmf() */
557
558 /*---------------------------------------------------------------------------*/
559
560 int
unur_distr_discr_set_cdf(struct unur_distr * distr,UNUR_FUNCT_DISCR * cdf)561 unur_distr_discr_set_cdf( struct unur_distr *distr, UNUR_FUNCT_DISCR *cdf )
562 /*----------------------------------------------------------------------*/
563 /* set CDF of distribution */
564 /* */
565 /* parameters: */
566 /* distr ... pointer to distribution object */
567 /* cdf ... pointer to CDF */
568 /* */
569 /* return: */
570 /* UNUR_SUCCESS ... on success */
571 /* error code ... on error */
572 /*----------------------------------------------------------------------*/
573 {
574 /* check arguments */
575 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
576 _unur_check_NULL( distr->name, cdf, UNUR_ERR_NULL );
577 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
578
579 /* it is not possible to set both a CDF and a PV */
580 if (DISTR.pv != NULL) {
581 _unur_warning(distr->name,UNUR_ERR_DISTR_SET,"delete exisiting PV");
582 free(DISTR.pv); DISTR.n_pv = 0;
583 }
584
585 /* we do not allow overwriting a CDF */
586 if (DISTR.cdf != NULL) {
587 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of CDF not allowed");
588 return UNUR_ERR_DISTR_SET;
589 }
590
591 /* changelog */
592 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
593 /* derived parameters like mode, sum, etc. might be wrong now! */
594
595 DISTR.cdf = cdf;
596 return UNUR_SUCCESS;
597 } /* end of unur_distr_discr_set_cdf() */
598
599 /*---------------------------------------------------------------------------*/
600
601 int
unur_distr_discr_set_invcdf(struct unur_distr * distr,UNUR_IFUNCT_DISCR * invcdf)602 unur_distr_discr_set_invcdf( struct unur_distr *distr, UNUR_IFUNCT_DISCR *invcdf )
603 /*----------------------------------------------------------------------*/
604 /* set inverse CDF of distribution */
605 /* */
606 /* parameters: */
607 /* distr ... pointer to distribution object */
608 /* invcdf ... pointer to inverse CDF */
609 /* */
610 /* return: */
611 /* UNUR_SUCCESS ... on success */
612 /* error code ... on error */
613 /*----------------------------------------------------------------------*/
614 {
615 /* check arguments */
616 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
617 _unur_check_NULL( distr->name, invcdf,UNUR_ERR_NULL );
618 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
619
620 /* we do not allow overwriting an inverse cdf */
621 if (DISTR.invcdf != NULL) {
622 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of inverse CDF not allowed");
623 return UNUR_ERR_DISTR_SET;
624 }
625
626 /* for derived distributions (e.g. order statistics) not possible */
627 if (distr->base) return UNUR_ERR_DISTR_INVALID;
628
629 /* changelog */
630 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
631 /* derived parameters like mode, area, etc. might be wrong now! */
632
633 DISTR.invcdf = invcdf;
634 return UNUR_SUCCESS;
635 } /* end of unur_distr_discr_set_invcdf() */
636
637 /*---------------------------------------------------------------------------*/
638
639 UNUR_FUNCT_DISCR *
unur_distr_discr_get_pmf(const struct unur_distr * distr)640 unur_distr_discr_get_pmf( const struct unur_distr *distr )
641 /*----------------------------------------------------------------------*/
642 /* get pointer to PMF of distribution */
643 /* */
644 /* parameters: */
645 /* distr ... pointer to distribution object */
646 /* */
647 /* return: */
648 /* pointer to PMF */
649 /*----------------------------------------------------------------------*/
650 {
651 /* check arguments */
652 _unur_check_NULL( NULL, distr, NULL );
653 _unur_check_distr_object( distr, DISCR, NULL );
654
655 return DISTR.pmf;
656 } /* end of unur_distr_discr_get_pmf() */
657
658 /*---------------------------------------------------------------------------*/
659
660 UNUR_FUNCT_DISCR *
unur_distr_discr_get_cdf(const struct unur_distr * distr)661 unur_distr_discr_get_cdf( const struct unur_distr *distr )
662 /*----------------------------------------------------------------------*/
663 /* get pointer to CDF of distribution */
664 /* */
665 /* parameters: */
666 /* distr ... pointer to distribution object */
667 /* */
668 /* return: */
669 /* pointer to CDF */
670 /*----------------------------------------------------------------------*/
671 {
672 /* check arguments */
673 _unur_check_NULL( NULL, distr, NULL );
674 _unur_check_distr_object( distr, DISCR, NULL );
675
676 return DISTR.cdf;
677 } /* end of unur_distr_discr_get_cdf() */
678
679 /*---------------------------------------------------------------------------*/
680
681 UNUR_IFUNCT_DISCR *
unur_distr_discr_get_invcdf(const struct unur_distr * distr)682 unur_distr_discr_get_invcdf( const struct unur_distr *distr )
683 /*----------------------------------------------------------------------*/
684 /* get pointer to inverse CDF of distribution */
685 /* */
686 /* parameters: */
687 /* distr ... pointer to distribution object */
688 /* */
689 /* return: */
690 /* pointer to inverse CDF */
691 /*----------------------------------------------------------------------*/
692 {
693 /* check arguments */
694 _unur_check_NULL( NULL, distr, NULL );
695 _unur_check_distr_object( distr, DISCR, NULL );
696
697 return DISTR.invcdf;
698 } /* end of unur_distr_discr_get_invcdf() */
699
700 /*---------------------------------------------------------------------------*/
701
702 double
unur_distr_discr_eval_pmf(int k,const struct unur_distr * distr)703 unur_distr_discr_eval_pmf( int k, const struct unur_distr *distr )
704 /*----------------------------------------------------------------------*/
705 /* evaluate PMF of distribution at k */
706 /* */
707 /* parameters: */
708 /* k ... argument for pmf */
709 /* distr ... pointer to distribution object */
710 /* */
711 /* return: */
712 /* pmf(k) */
713 /*----------------------------------------------------------------------*/
714 {
715 /* check arguments */
716 _unur_check_NULL( NULL, distr, INFINITY );
717 _unur_check_distr_object( distr, DISCR, INFINITY );
718
719 if (DISTR.pmf == NULL) {
720 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
721 return INFINITY;
722 }
723
724 return _unur_discr_PMF(k,distr);
725 } /* end of unur_distr_discr_eval_pmf() */
726
727 /*---------------------------------------------------------------------------*/
728
729 double
unur_distr_discr_eval_cdf(int k,const struct unur_distr * distr)730 unur_distr_discr_eval_cdf( int k, const struct unur_distr *distr )
731 /*----------------------------------------------------------------------*/
732 /* evaluate CDF of distribution at k */
733 /* */
734 /* parameters: */
735 /* k ... argument for CDF */
736 /* distr ... pointer to distribution object */
737 /* */
738 /* return: */
739 /* CDF(k) */
740 /*----------------------------------------------------------------------*/
741 {
742 /* check arguments */
743 _unur_check_NULL( NULL, distr, INFINITY );
744 _unur_check_distr_object( distr, DISCR, INFINITY );
745
746 if (DISTR.cdf == NULL) {
747 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
748 return INFINITY;
749 }
750
751 return _unur_discr_CDF(k,distr);
752 } /* end of unur_distr_discr_eval_cdf() */
753
754 /*---------------------------------------------------------------------------*/
755
756 int
unur_distr_discr_eval_invcdf(double u,const struct unur_distr * distr)757 unur_distr_discr_eval_invcdf( double u, const struct unur_distr *distr )
758 /*----------------------------------------------------------------------*/
759 /* evaluate inverse CDF of distribution at u */
760 /* */
761 /* parameters: */
762 /* u ... argument for inverse CDF */
763 /* distr ... pointer to distribution object */
764 /* */
765 /* return: */
766 /* invcdf(u) */
767 /* INT_MAX ... on error */
768 /*----------------------------------------------------------------------*/
769 {
770 /* check arguments */
771 _unur_check_NULL( NULL, distr, INT_MAX );
772 _unur_check_distr_object( distr, DISCR, INT_MAX );
773
774 if (DISTR.invcdf == NULL) {
775 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
776 return INT_MAX;
777 }
778
779 if (u<=0.)
780 return DISTR.domain[0];
781 if (u>=1.)
782 return DISTR.domain[1];
783 else
784 return _unur_discr_invCDF(u,distr);
785
786 } /* end of unur_distr_discr_eval_invcdf() */
787
788 /*---------------------------------------------------------------------------*/
789
790 int
unur_distr_discr_set_pmfstr(struct unur_distr * distr,const char * pmfstr)791 unur_distr_discr_set_pmfstr( struct unur_distr *distr, const char *pmfstr )
792 /*----------------------------------------------------------------------*/
793 /* set PMF of distribution via a string interface */
794 /* */
795 /* parameters: */
796 /* distr ... pointer to distribution object */
797 /* pmfstr ... string that describes function term of PMF */
798 /* */
799 /* return: */
800 /* UNUR_SUCCESS ... on success */
801 /* error code ... on error */
802 /*----------------------------------------------------------------------*/
803 {
804 /* check arguments */
805 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
806 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
807 _unur_check_NULL( NULL, pmfstr, UNUR_ERR_NULL );
808
809 /* it is not possible to set a PMF when a PV is given. */
810 if (DISTR.pv != NULL) {
811 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"PV given, cannot set PMF");
812 return UNUR_ERR_DISTR_SET;
813 }
814
815 /* we do not allow overwriting a PMF */
816 if (DISTR.pmf != NULL) {
817 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of PMF not allowed");
818 return UNUR_ERR_DISTR_SET;
819 }
820
821 /* for derived distributions (e.g. order statistics) not possible */
822 if (distr->base) return UNUR_ERR_DISTR_DATA;
823
824 /* changelog */
825 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
826 /* derived parameters like mode, area, etc. might be wrong now! */
827
828 /* parse PMF string */
829 if ( (DISTR.pmftree = _unur_fstr2tree(pmfstr)) == NULL ) {
830 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Syntax error in function string");
831 return UNUR_ERR_DISTR_SET;
832 }
833 DISTR.pmf = _unur_distr_discr_eval_pmf_tree;
834
835 return UNUR_SUCCESS;
836 } /* end of unur_distr_discr_set_pmfstr() */
837
838 /*---------------------------------------------------------------------------*/
839
840 int
unur_distr_discr_set_cdfstr(struct unur_distr * distr,const char * cdfstr)841 unur_distr_discr_set_cdfstr( struct unur_distr *distr, const char *cdfstr )
842 /*----------------------------------------------------------------------*/
843 /* set CDF of distribution via a string interface */
844 /* */
845 /* parameters: */
846 /* distr ... pointer to distribution object */
847 /* cdfstr ... string that describes function term of CDF */
848 /* */
849 /* return: */
850 /* UNUR_SUCCESS ... on success */
851 /* error code ... on error */
852 /*----------------------------------------------------------------------*/
853 {
854 /* check arguments */
855 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
856 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
857 _unur_check_NULL( NULL, cdfstr, UNUR_ERR_NULL );
858
859 /* we do not allow overwriting a CDF */
860 if (DISTR.cdf != NULL) {
861 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Overwriting of CDF not allowed");
862 return UNUR_ERR_DISTR_SET;
863 }
864
865 /* for derived distributions (e.g. order statistics) not possible */
866 if (distr->base) return UNUR_ERR_DISTR_DATA;
867
868 /* changelog */
869 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
870 /* derived parameters like mode, area, etc. might be wrong now! */
871
872 /* parse string */
873 if ( (DISTR.cdftree = _unur_fstr2tree(cdfstr)) == NULL ) {
874 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"Syntax error in function string");
875 return UNUR_ERR_DISTR_SET;
876 }
877
878 /* set evaluation function */
879 DISTR.cdf = _unur_distr_discr_eval_cdf_tree;
880
881 return UNUR_SUCCESS;
882 } /* end of unur_distr_discr_set_cdfstr() */
883
884 /*---------------------------------------------------------------------------*/
885
886 double
_unur_distr_discr_eval_pmf_tree(int k,const struct unur_distr * distr)887 _unur_distr_discr_eval_pmf_tree( int k, const struct unur_distr *distr )
888 /*----------------------------------------------------------------------*/
889 /* evaluate function tree for PMF. */
890 /* */
891 /* parameters: */
892 /* k ... argument for PMF */
893 /* distr ... pointer to distribution object */
894 /* */
895 /* return: */
896 /* PMF at k */
897 /*----------------------------------------------------------------------*/
898 {
899 /* check arguments */
900 _unur_check_NULL( NULL, distr, INFINITY );
901 _unur_check_distr_object( distr, DISCR, INFINITY );
902
903 return ((DISTR.pmftree) ? _unur_fstr_eval_tree(DISTR.pmftree,(double)k) : 0.);
904 } /* end of _unur_distr_discr_eval_pmf_tree() */
905
906 /*---------------------------------------------------------------------------*/
907
908 double
_unur_distr_discr_eval_cdf_tree(int k,const struct unur_distr * distr)909 _unur_distr_discr_eval_cdf_tree( int k, const struct unur_distr *distr )
910 /*----------------------------------------------------------------------*/
911 /* evaluate function tree for CDF. */
912 /* */
913 /* parameters: */
914 /* k ... argument for CDF */
915 /* distr ... pointer to distribution object */
916 /* */
917 /* return: */
918 /* CDF at k */
919 /*----------------------------------------------------------------------*/
920 {
921 /* check arguments */
922 _unur_check_NULL( NULL, distr, INFINITY );
923 _unur_check_distr_object( distr, DISCR, INFINITY );
924
925 return ((DISTR.cdftree) ? _unur_fstr_eval_tree(DISTR.cdftree,(double)k) : 0.);
926 } /* end of _unur_distr_discr_eval_cdf_tree() */
927
928 /*---------------------------------------------------------------------------*/
929
930 char *
unur_distr_discr_get_pmfstr(const struct unur_distr * distr)931 unur_distr_discr_get_pmfstr( const struct unur_distr *distr )
932 /*----------------------------------------------------------------------*/
933 /* get PMF string that is given via the string interface */
934 /* */
935 /* parameters: */
936 /* distr ... pointer to distribution object */
937 /* */
938 /* return: */
939 /* pointer to resulting string. */
940 /* */
941 /* comment: */
942 /* This string should be freed when it is not used any more. */
943 /*----------------------------------------------------------------------*/
944 {
945 /* check arguments */
946 _unur_check_NULL( NULL, distr, NULL );
947 _unur_check_distr_object( distr, DISCR, NULL );
948 _unur_check_NULL( NULL, DISTR.pmftree, NULL );
949
950 /* make and return string */
951 return _unur_fstr_tree2string(DISTR.pmftree,"x","PMF",TRUE);
952 } /* end of unur_distr_discr_get_pmfstr() */
953
954 /*---------------------------------------------------------------------------*/
955
956 char *
unur_distr_discr_get_cdfstr(const struct unur_distr * distr)957 unur_distr_discr_get_cdfstr( const struct unur_distr *distr )
958 /*----------------------------------------------------------------------*/
959 /* get CDF string that is given via the string interface */
960 /* */
961 /* parameters: */
962 /* distr ... pointer to distribution object */
963 /* */
964 /* return: */
965 /* pointer to resulting string. */
966 /* */
967 /* comment: */
968 /* This string should be freed when it is not used any more. */
969 /*----------------------------------------------------------------------*/
970 {
971 /* check arguments */
972 _unur_check_NULL( NULL, distr, NULL );
973 _unur_check_distr_object( distr, DISCR, NULL );
974 _unur_check_NULL( NULL, DISTR.cdftree, NULL );
975
976 /* make and return string */
977 return _unur_fstr_tree2string(DISTR.cdftree,"x","CDF",TRUE);
978 } /* end of unur_distr_discr_get_cdfstr() */
979
980 /*---------------------------------------------------------------------------*/
981
982 int
unur_distr_discr_set_pmfparams(struct unur_distr * distr,const double * params,int n_params)983 unur_distr_discr_set_pmfparams( struct unur_distr *distr, const double *params, int n_params )
984 /*----------------------------------------------------------------------*/
985 /* set array of parameters for distribution */
986 /* */
987 /* parameters: */
988 /* distr ... pointer to distribution object */
989 /* params ... list of arguments */
990 /* n_params ... number of arguments */
991 /* */
992 /* return: */
993 /* UNUR_SUCCESS ... on success */
994 /* error code ... on error */
995 /*----------------------------------------------------------------------*/
996 {
997 /* check arguments */
998 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
999 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
1000 if (n_params>0) _unur_check_NULL(distr->name, params, UNUR_ERR_NULL);
1001
1002 /* first check number of new parameter for the distribution */
1003 if (n_params < 0 || n_params > UNUR_DISTR_MAXPARAMS ) {
1004 _unur_error(NULL,UNUR_ERR_DISTR_NPARAMS,"");
1005 return UNUR_ERR_DISTR_NPARAMS;
1006 }
1007
1008 /* changelog */
1009 distr->set &= ~UNUR_DISTR_SET_MASK_DERIVED;
1010 /* derived parameters like mode, area, etc. might be wrong now! */
1011
1012 /* even if the set routine fails, the derived parameters are
1013 marked as unknown. but this is o.k. since in this case something
1014 has been wrong. */
1015
1016 /* use special routine for setting parameters
1017 (if there is one) */
1018
1019 if (DISTR.set_params)
1020 return (DISTR.set_params(distr,params,n_params));
1021
1022 /* otherwise simply copy parameters */
1023
1024 DISTR.n_params = n_params;
1025 if (n_params) memcpy( DISTR.params, params, n_params*sizeof(double) );
1026
1027 /* o.k. */
1028 return UNUR_SUCCESS;
1029 } /* end of unur_distr_discr_set_pmfparams() */
1030
1031 /*---------------------------------------------------------------------------*/
1032
1033 int
unur_distr_discr_get_pmfparams(const struct unur_distr * distr,const double ** params)1034 unur_distr_discr_get_pmfparams( const struct unur_distr *distr, const double **params )
1035 /*----------------------------------------------------------------------*/
1036 /* get number of pmf parameters and sets pointer to array params[] of */
1037 /* parameters */
1038 /* */
1039 /* parameters: */
1040 /* distr ... pointer to distribution object */
1041 /* params ... pointer to list of arguments */
1042 /* */
1043 /* return: */
1044 /* number of pmf parameters */
1045 /* */
1046 /* error: */
1047 /* return 0 */
1048 /*----------------------------------------------------------------------*/
1049 {
1050 /* check arguments */
1051 _unur_check_NULL( NULL, distr, 0 );
1052 _unur_check_distr_object( distr, DISCR, 0 );
1053
1054 *params = (DISTR.n_params) ? DISTR.params : NULL;
1055 return DISTR.n_params;
1056
1057 } /* end of unur_distr_discr_get_pmfparams() */
1058
1059 /*---------------------------------------------------------------------------*/
1060
1061 int
unur_distr_discr_set_domain(struct unur_distr * distr,int left,int right)1062 unur_distr_discr_set_domain( struct unur_distr *distr, int left, int right )
1063 /*----------------------------------------------------------------------*/
1064 /* set the left and right borders of the domain of the distribution */
1065 /* */
1066 /* parameters: */
1067 /* distr ... pointer to distribution object */
1068 /* left ... left boundary point */
1069 /* right ... right boundary point */
1070 /* */
1071 /* return: */
1072 /* UNUR_SUCCESS ... on success */
1073 /* error code ... on error */
1074 /* */
1075 /* comment: */
1076 /* INT_MIN and INT_MAX are interpreted as (minus) infinity */
1077 /*----------------------------------------------------------------------*/
1078 {
1079 /* check arguments */
1080 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1081 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
1082
1083 /* check new parameter for distribution */
1084 if (left >= right) {
1085 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"domain, left >= right");
1086 return UNUR_ERR_DISTR_SET;
1087 }
1088
1089 /* store data */
1090 DISTR.trunc[0] = DISTR.domain[0] = left;
1091 DISTR.trunc[1] = DISTR.domain[1] = (DISTR.pv == NULL) ? right : left+DISTR.n_pv-1;
1092
1093 /* changelog */
1094 distr->set |= UNUR_DISTR_SET_DOMAIN;
1095
1096 /* if distr is an object for a standard distribution, this */
1097 /* not the original domain of it. (not a "standard domain") */
1098 /* However, since we have changed the domain, we assume */
1099 /* that this is not a truncated distribution. */
1100 /* At last we have to mark all derived parameters as unknown */
1101 distr->set &= ~(UNUR_DISTR_SET_STDDOMAIN |
1102 UNUR_DISTR_SET_TRUNCATED |
1103 UNUR_DISTR_SET_MASK_DERIVED );
1104
1105 /* o.k. */
1106 return UNUR_SUCCESS;
1107
1108 } /* end of unur_distr_discr_set_domain() */
1109
1110 /*---------------------------------------------------------------------------*/
1111
1112 int
unur_distr_discr_get_domain(const struct unur_distr * distr,int * left,int * right)1113 unur_distr_discr_get_domain( const struct unur_distr *distr, int *left, int *right )
1114 /*----------------------------------------------------------------------*/
1115 /* set the left and right borders of the domain of the distribution */
1116 /* */
1117 /* parameters: */
1118 /* distr ... pointer to distribution object */
1119 /* left ... left boundary point */
1120 /* right ... right boundary point */
1121 /* */
1122 /* return: */
1123 /* UNUR_SUCCESS ... on success */
1124 /* error code ... on error */
1125 /* */
1126 /* comment: */
1127 /* INT_MIN and INT_MAX are interpreted as (minus) infinity */
1128 /* if no boundaries have been set [INT_MIN, INT_MAX] is returned. */
1129 /*----------------------------------------------------------------------*/
1130 {
1131 /* in case of error the boundaries are set to +/- INFINITY */
1132 *left = INT_MIN;
1133 *right = INT_MAX;
1134
1135 /* check arguments */
1136 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1137 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
1138
1139 /* o.k. */
1140 *left = DISTR.domain[0];
1141 *right = DISTR.domain[1];
1142
1143 return UNUR_SUCCESS;
1144 } /* end of unur_distr_discr_get_domain() */
1145
1146 /*---------------------------------------------------------------------------*/
1147
1148 int
unur_distr_discr_set_mode(struct unur_distr * distr,int mode)1149 unur_distr_discr_set_mode( struct unur_distr *distr, int mode )
1150 /*----------------------------------------------------------------------*/
1151 /* set mode of distribution */
1152 /* */
1153 /* parameters: */
1154 /* distr ... pointer to distribution object */
1155 /* mode ... mode of PMF */
1156 /* */
1157 /* return: */
1158 /* UNUR_SUCCESS ... on success */
1159 /* error code ... on error */
1160 /*----------------------------------------------------------------------*/
1161 {
1162 /* check arguments */
1163 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1164 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
1165
1166 DISTR.mode = mode;
1167
1168 /* changelog */
1169 distr->set |= UNUR_DISTR_SET_MODE;
1170
1171 /* o.k. */
1172 return UNUR_SUCCESS;
1173 } /* end of unur_distr_discr_set_mode() */
1174
1175 /*---------------------------------------------------------------------------*/
1176
1177 int
unur_distr_discr_upd_mode(struct unur_distr * distr)1178 unur_distr_discr_upd_mode( struct unur_distr *distr )
1179 /*----------------------------------------------------------------------*/
1180 /* (re-) compute mode of distribution (if possible) */
1181 /* */
1182 /* parameters: */
1183 /* distr ... pointer to distribution object */
1184 /* */
1185 /* return: */
1186 /* UNUR_SUCCESS ... on success */
1187 /* error code ... on error */
1188 /*----------------------------------------------------------------------*/
1189 {
1190 /* check arguments */
1191 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1192 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
1193
1194 if (DISTR.upd_mode == NULL) {
1195 /* no function to compute mode available */
1196 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
1197 return UNUR_ERR_DISTR_DATA;
1198 }
1199
1200 /* compute mode */
1201 if ((DISTR.upd_mode)(distr)==UNUR_SUCCESS) {
1202 /* changelog */
1203 distr->set |= UNUR_DISTR_SET_MODE;
1204 return UNUR_SUCCESS;
1205 }
1206 else {
1207 /* computing of mode failed */
1208 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"");
1209 return UNUR_ERR_DISTR_DATA;
1210 }
1211
1212 } /* end of unur_distr_discr_upd_mode() */
1213
1214 /*---------------------------------------------------------------------------*/
1215
1216 int
unur_distr_discr_get_mode(struct unur_distr * distr)1217 unur_distr_discr_get_mode( struct unur_distr *distr )
1218 /*----------------------------------------------------------------------*/
1219 /* get mode of distribution */
1220 /* */
1221 /* parameters: */
1222 /* distr ... pointer to distribution object */
1223 /* */
1224 /* return: */
1225 /* mode of distribution */
1226 /*----------------------------------------------------------------------*/
1227 {
1228 /* check arguments */
1229 _unur_check_NULL( NULL, distr, INT_MAX );
1230 _unur_check_distr_object( distr, DISCR, INT_MAX );
1231
1232 /* mode known ? */
1233 if ( !(distr->set & UNUR_DISTR_SET_MODE) ) {
1234 /* try to compute mode */
1235 if (DISTR.upd_mode == NULL) {
1236 /* no function to compute mode available */
1237 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"mode");
1238 return INT_MAX;
1239 }
1240 else {
1241 /* compute mode */
1242 if (unur_distr_discr_upd_mode(distr)!=UNUR_SUCCESS) {
1243 /* finding mode not successfully */
1244 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"mode");
1245 return INT_MAX;
1246 }
1247 }
1248 }
1249
1250 return DISTR.mode;
1251
1252 } /* end of unur_distr_discr_get_mode() */
1253
1254 /*---------------------------------------------------------------------------*/
1255
1256 int
unur_distr_discr_set_pmfsum(struct unur_distr * distr,double sum)1257 unur_distr_discr_set_pmfsum( struct unur_distr *distr, double sum )
1258 /*----------------------------------------------------------------------*/
1259 /* set sum over PMF */
1260 /* */
1261 /* parameters: */
1262 /* distr ... pointer to distribution object */
1263 /* sum ... sum over PMF */
1264 /* */
1265 /* return: */
1266 /* UNUR_SUCCESS ... on success */
1267 /* error code ... on error */
1268 /*----------------------------------------------------------------------*/
1269 {
1270 /* check arguments */
1271 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1272 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
1273
1274 /* check new parameter for distribution */
1275 if (sum <= 0.) {
1276 _unur_error(distr->name,UNUR_ERR_DISTR_SET,"pmf sum <= 0");
1277 return UNUR_ERR_DISTR_SET;
1278 }
1279
1280 DISTR.sum = sum;
1281
1282 /* changelog */
1283 distr->set |= UNUR_DISTR_SET_PMFSUM;
1284
1285 /* o.k. */
1286 return UNUR_SUCCESS;
1287
1288 } /* end of unur_distr_discr_set_pmfsum() */
1289
1290 /*---------------------------------------------------------------------------*/
1291
1292 int
unur_distr_discr_upd_pmfsum(struct unur_distr * distr)1293 unur_distr_discr_upd_pmfsum( struct unur_distr *distr )
1294 /*----------------------------------------------------------------------*/
1295 /* (re-) compute sum over PMF of distribution (if possible) */
1296 /* */
1297 /* parameters: */
1298 /* distr ... pointer to distribution object */
1299 /* */
1300 /* return: */
1301 /* UNUR_SUCCESS ... on success */
1302 /* error code ... on error */
1303 /*----------------------------------------------------------------------*/
1304 {
1305 double sum = 0.;
1306 int k, left, right, length;
1307
1308 /* check arguments */
1309 _unur_check_NULL( NULL, distr, UNUR_ERR_NULL );
1310 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_SET );
1311
1312 /* changelog */
1313 distr->set |= UNUR_DISTR_SET_PMFSUM;
1314
1315 if (DISTR.upd_sum != NULL) {
1316 /* try function given by distribution */
1317 if ((DISTR.upd_sum)(distr)==UNUR_SUCCESS)
1318 return UNUR_SUCCESS;
1319 }
1320
1321 /* no function to compute sum available */
1322 left = DISTR.domain[0];
1323 right = DISTR.domain[1];
1324 length = right - left;
1325 /* remark: length < 0 if right-left overflows */
1326
1327
1328 if (DISTR.cdf != NULL) {
1329 /* use CDF */
1330 if (left > INT_MIN) left -= 1;
1331 DISTR.sum = _unur_discr_CDF(right,distr) - _unur_discr_CDF(left,distr);
1332 return UNUR_SUCCESS;
1333 }
1334
1335 if (DISTR.pv != NULL) {
1336 for (k = 0; k<= length; k++)
1337 /* use probability vector */
1338 sum += DISTR.pv[k];
1339 DISTR.sum = sum;
1340 return UNUR_SUCCESS;
1341 }
1342
1343 if (DISTR.pmf != NULL && length > 0 && length <= MAX_PMF_DOMAIN_FOR_UPD_PMFSUM) {
1344 /* use PMF */
1345 for (k = left; k<= right; k++)
1346 sum += _unur_discr_PMF(k,distr);
1347 DISTR.sum = sum;
1348 return UNUR_SUCCESS;
1349 }
1350
1351 /* error: data missing */
1352 distr->set &= ~UNUR_DISTR_SET_PMFSUM;
1353 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,"Cannot compute sum");
1354 return UNUR_ERR_DISTR_DATA;
1355
1356 } /* end of unur_distr_discr_upd_pmfsum() */
1357
1358 /*---------------------------------------------------------------------------*/
1359
1360 double
unur_distr_discr_get_pmfsum(struct unur_distr * distr)1361 unur_distr_discr_get_pmfsum( struct unur_distr *distr )
1362 /*----------------------------------------------------------------------*/
1363 /* get sum over PMF of distribution */
1364 /* */
1365 /* parameters: */
1366 /* distr ... pointer to distribution object */
1367 /* */
1368 /* return: */
1369 /* sum over PMF of distribution */
1370 /*----------------------------------------------------------------------*/
1371 {
1372 /* check arguments */
1373 _unur_check_NULL( NULL, distr, INFINITY );
1374 _unur_check_distr_object( distr, DISCR, INFINITY );
1375
1376 /* sum known ? */
1377 if ( !(distr->set & UNUR_DISTR_SET_PMFSUM) ) {
1378 /* try to compute sum */
1379 if ( unur_distr_discr_upd_pmfsum(distr) != UNUR_SUCCESS ) {
1380 _unur_error(distr->name,UNUR_ERR_DISTR_GET,"sum");
1381 return INFINITY;
1382 }
1383 }
1384
1385 return DISTR.sum;
1386
1387 } /* end of unur_distr_discr_get_pmfsum() */
1388
1389 /*****************************************************************************/
1390
1391 /*---------------------------------------------------------------------------*/
1392 #ifdef UNUR_ENABLE_LOGGING
1393 /*---------------------------------------------------------------------------*/
1394
1395 void
_unur_distr_discr_debug(const struct unur_distr * distr,const char * genid,unsigned printvector)1396 _unur_distr_discr_debug( const struct unur_distr *distr, const char *genid, unsigned printvector )
1397 /*----------------------------------------------------------------------*/
1398 /* write info about distribution into LOG file */
1399 /* */
1400 /* parameters: */
1401 /* distr ... pointer to distribution object */
1402 /* genid ... pointer to generator id */
1403 /* printvector ... whether the probability vector (if given) */
1404 /*----------------------------------------------------------------------*/
1405 {
1406 FILE *LOG;
1407 int i;
1408
1409 /* check arguments */
1410 CHECK_NULL(distr,RETURN_VOID);
1411 COOKIE_CHECK(distr,CK_DISTR_DISCR,RETURN_VOID);
1412
1413 LOG = unur_get_stream();
1414
1415 fprintf(LOG,"%s: distribution:\n",genid);
1416 fprintf(LOG,"%s:\ttype = discrete univariate distribution\n",genid);
1417 fprintf(LOG,"%s:\tname = %s\n",genid,distr->name);
1418
1419 if ( DISTR.pmf ) {
1420 /* have probability mass function */
1421 fprintf(LOG,"%s:\tPMF with %d argument(s)\n",genid,DISTR.n_params);
1422 for( i=0; i<DISTR.n_params; i++ )
1423 fprintf(LOG,"%s:\t\tparam[%d] = %g\n",genid,i,DISTR.params[i]);
1424 }
1425
1426 if (DISTR.n_pv>0) {
1427 /* have probability vector */
1428 fprintf(LOG,"%s:\tprobability vector of length %d",genid,DISTR.n_pv);
1429 if (printvector) {
1430 for (i=0; i<DISTR.n_pv; i++) {
1431 if (i%10 == 0)
1432 fprintf(LOG,"\n%s:\t",genid);
1433 fprintf(LOG," %.5f",DISTR.pv[i]);
1434 }
1435 }
1436 fprintf(LOG,"\n%s:\n",genid);
1437 }
1438
1439 /* domain */
1440 if ( DISTR.pmf ) {
1441 /* have probability mass function */
1442 fprintf(LOG,"%s:\tdomain for pmf = (%d, %d)",genid,DISTR.domain[0],DISTR.domain[1]);
1443 _unur_print_if_default(distr,UNUR_DISTR_SET_DOMAIN);
1444 fprintf(LOG,"\n%s:\n",genid);
1445 }
1446
1447 if (DISTR.n_pv>0) {
1448 /* have probability vector */
1449 fprintf(LOG,"%s:\tdomain for pv = (%d, %d)",genid,DISTR.domain[0],DISTR.domain[0]-1+DISTR.n_pv);
1450 _unur_print_if_default(distr,UNUR_DISTR_SET_DOMAIN);
1451 fprintf(LOG,"\n%s:\n",genid);
1452 }
1453
1454 if (distr->set & UNUR_DISTR_SET_MODE)
1455 fprintf(LOG,"%s:\tmode = %d\n",genid,DISTR.mode);
1456 else
1457 fprintf(LOG,"%s:\tmode unknown\n",genid);
1458
1459 if (distr->set & UNUR_DISTR_SET_PMFSUM)
1460 fprintf(LOG,"%s:\tsum over PMF = %g\n",genid,DISTR.sum);
1461 else
1462 fprintf(LOG,"%s:\tsum over PMF unknown\n",genid);
1463 fprintf(LOG,"%s:\n",genid);
1464
1465 } /* end of _unur_distr_discr_debug() */
1466
1467 /*---------------------------------------------------------------------------*/
1468 #endif /* end UNUR_ENABLE_LOGGING */
1469 /*---------------------------------------------------------------------------*/
1470
1471 /*****************************************************************************/
1472
1473 /*---------------------------------------------------------------------------*/
1474
1475 int
_unur_distr_discr_find_mode(struct unur_distr * distr)1476 _unur_distr_discr_find_mode(struct unur_distr *distr )
1477 /*----------------------------------------------------------------------*/
1478 /* find mode of a probability vector by bisection. */
1479 /* */
1480 /* Assumptions: PMF is unimodal, no "inflexion" ("saddle") points. */
1481 /* */
1482 /* return: */
1483 /* UNUR_SUCCESS ... on success */
1484 /* error code ... on error */
1485 /*----------------------------------------------------------------------*/
1486 {
1487 #define N_TRIALS (100)
1488
1489 int x[3]; /* bracket for mode x[0] < x[1] < x[2] */
1490 double fx[3]; /* PMF values at bracket points */
1491 int xnew; /* new point in iteration */
1492 double fxnew; /* PMF value at new point */
1493 int step; /* auxiliary variable for searching point */
1494 int this, other; /* side of bracket under consideration */
1495 int cutthis; /* which side of bracket should be cut */
1496
1497 const double r = (sqrt(5.)-1.)/2.; /* sectio aurea */
1498
1499 /* check arguments */
1500 CHECK_NULL( distr, UNUR_ERR_NULL );
1501 _unur_check_distr_object( distr, DISCR, UNUR_ERR_DISTR_INVALID );
1502
1503 /* left and right boundary point of bracket */
1504 x[0] = DISTR.domain[0];
1505 x[2] = DISTR.domain[1];
1506 fx[0] = unur_distr_discr_eval_pv(x[0], distr);
1507 fx[2] = unur_distr_discr_eval_pv(x[2], distr);
1508
1509 /* we assume for our algorithm that there are at least three points */
1510 if (x[2] <= x[0] + 1) {
1511 /* we have one or two points only */
1512 DISTR.mode = (fx[0] <= fx[2]) ? x[2] : x[0];
1513 distr->set |= UNUR_DISTR_SET_MODE | UNUR_DISTR_SET_MODE_APPROX ;
1514 return UNUR_SUCCESS;
1515 }
1516
1517 /* --- middle point of bracket ------------------------------------------- */
1518 /* we have to find a middle point where the PMF is non-zero */
1519
1520 /* first trial: mean of boundary points */
1521 x[1] = (x[0]/2) + (x[2]/2);
1522 if (x[1]<=x[0]) x[1]++;
1523 if (x[1]>=x[2]) x[1]--;
1524 fx[1] = unur_distr_discr_eval_pv(x[1], distr);
1525
1526 /* second trial: start search from left boundary */
1527 if ( !(fx[1]>0.)) {
1528 xnew = (DISTR.domain[0]!=INT_MIN) ? DISTR.domain[0] : 0;
1529 for (step = 1; step < N_TRIALS; step++) {
1530 xnew += step;
1531 if (xnew >= DISTR.domain[1]) break;
1532 if ((fxnew = unur_distr_discr_eval_pv(xnew,distr)) > 0.) {
1533 x[1] = xnew; fx[1] = fxnew; break;
1534 }
1535 }
1536 }
1537
1538 /* third trial: start search from 0 */
1539 if ( !(fx[1]>0.) && DISTR.domain[0]!=0) {
1540 xnew = 0;
1541 for (step = 1; step < N_TRIALS; step++) {
1542 xnew += step;
1543 if (xnew >= DISTR.domain[1]) break;
1544 if ((fxnew = unur_distr_discr_eval_pv(xnew,distr)) > 0.) {
1545 x[1] = xnew; fx[1] = fxnew; break;
1546 }
1547 }
1548 }
1549
1550 /* forth trial: start search from right boundary */
1551 if ( !(fx[1]>0.) && DISTR.domain[1]!=INT_MAX) {
1552 xnew = DISTR.domain[1];
1553 for (step = 1; step < N_TRIALS; step++) {
1554 xnew -= step;
1555 if (xnew <= DISTR.domain[0]) break;
1556 if ((fxnew = unur_distr_discr_eval_pv(xnew,distr)) > 0.) {
1557 x[1] = xnew; fx[1] = fxnew; break;
1558 }
1559 }
1560 }
1561
1562 if ( !(fx[1]>0.)) {
1563 _unur_error(distr->name,UNUR_ERR_DISTR_DATA,
1564 "find_mode(): no positive entry in PV found");
1565 return UNUR_ERR_DISTR_DATA;
1566 }
1567 if (fx[1]<fx[0] && fx[1]<fx[2]) {
1568 _unur_error(distr->name,UNUR_ERR_DISTR_DATA, "find_mode(): PV not unimodal");
1569 return UNUR_ERR_DISTR_DATA;
1570 }
1571
1572 /* --- middle point of bracket found ------------------------------------- */
1573
1574 /* --- interate until maximum is found ----------------------------------- */
1575
1576 while (1) {
1577
1578 /* fprintf(stderr,"x = %d, %d, %d\n",x[0],x[1],x[2]); */
1579 /* fprintf(stderr,"fx = %g, %g, %g\n",fx[0],fx[1],fx[2]); */
1580
1581 /* terminate */
1582 if (x[0]+1 >= x[1] && x[1] >= x[2]-1) {
1583 DISTR.mode = (fx[0]>fx[2]) ? x[0] : x[2];
1584 if (fx[1]>DISTR.mode) DISTR.mode = x[1];
1585 distr->set |= UNUR_DISTR_SET_MODE | UNUR_DISTR_SET_MODE_APPROX ;
1586 return UNUR_SUCCESS;
1587 }
1588
1589 /* new point */
1590 xnew = (int) (r*x[0] + (1.-r)*x[2]);
1591 if (xnew == x[0]) ++xnew;
1592 if (xnew == x[2]) --xnew;
1593 if (xnew == x[1]) xnew += (x[1]-1==x[0]) ? 1 : -1;
1594
1595 /* side of bracket */
1596 if (xnew < x[1]) {
1597 this = 0; other = 2; } /* l.h.s. of bracket */
1598 else {
1599 this = 2; other = 0; } /* r.h.s. of bracket */
1600
1601 /* value at new point */
1602 fxnew = unur_distr_discr_eval_pv(xnew,distr);
1603 if ( fxnew < fx[0] && fxnew < fx[2] ) {
1604 _unur_error(distr->name,UNUR_ERR_DISTR_DATA, "find_mode(): PV not unimodal");
1605 return UNUR_ERR_DISTR_DATA;
1606 }
1607
1608 do {
1609
1610 if (!_unur_FP_same(fxnew,fx[1])) {
1611 cutthis = (fxnew > fx[1]) ? FALSE : TRUE;
1612 break;
1613 }
1614
1615 /* else: fxnew == fx[1] */
1616
1617 if (fx[this] > fx[1]) { cutthis = FALSE; break; }
1618 if (fx[other] > fx[1]) { cutthis = TRUE; break; }
1619
1620 /* else: fx[0] < fxnew && fx[1] == fxnew && fx[2] < fxnew */
1621
1622 for (step = 1; step < N_TRIALS && xnew >= x[0] && xnew <= x[2]; step++) {
1623 xnew += (this==0) ? -1 : 1;
1624 fxnew = unur_distr_discr_eval_pv(xnew,distr);
1625 if (_unur_FP_less(fxnew,fx[1])) {
1626 DISTR.mode = x[1];
1627 distr->set |= UNUR_DISTR_SET_MODE | UNUR_DISTR_SET_MODE_APPROX ;
1628 return UNUR_SUCCESS;
1629 }
1630 }
1631
1632 _unur_error(distr->name,UNUR_ERR_DISTR_DATA, "find_mode(): PV not unimodal");
1633 return UNUR_ERR_DISTR_DATA;
1634
1635 } while (0);
1636
1637 if (cutthis) {
1638 x[this] = xnew; fx[this] = fxnew;
1639 }
1640 else {
1641 x[other] = x[1]; fx[other] = fx[1];
1642 x[1] = xnew; fx[1] = fxnew;
1643 }
1644
1645 } /* --- end while(1) --- */
1646
1647 /* changelog */
1648 /* distr->set |= UNUR_DISTR_SET_MODE | UNUR_DISTR_SET_MODE_APPROX ; */
1649
1650 /* o.k. */
1651 /* return UNUR_SUCCESS; */
1652
1653 } /* end of _unur_distr_discr_find_mode() */
1654
1655 /*---------------------------------------------------------------------------*/
1656 #undef DISTR
1657 /*---------------------------------------------------------------------------*/
1658