1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: c_burr.c *
8 * *
9 * REFERENCES: *
10 * *
11 * [2] N.L. Johnson, S. Kotz and N. Balakrishnan *
12 * Continuous Univariate Distributions, *
13 * Volume 1, 2nd edition *
14 * John Wiley & Sons, Inc., New York, 1994 *
15 * *
16 *****************************************************************************
17 * *
18 * Burr family of distributions [2; ch.12, p.54] *
19 * *
20 * first parameter must be type of distribution (1 ... 12)! *
21 * *
22 *...........................................................................*
23 * *
24 * distr: Burr distribution Type I *
25 * first parameter is 1 *
26 * *
27 * cdf: F(x) = x *
28 * pdf: f(x) = 1 *
29 * domain: 0 <= x <= 1 *
30 * constant: 1 *
31 * *
32 * parameters: none *
33 * *
34 *...........................................................................*
35 * *
36 * distr: Burr distribution Type II *
37 * first parameter is 2 *
38 * *
39 * cdf: F(x) = (exp(-x) + 1)^(-k) *
40 * pdf: f(x) = k * exp(-x) * (exp(-x)+1)^(-k-1) *
41 * domain: -infinity < x < infinity *
42 * constant: 1 *
43 * *
44 * parameters: 1 *
45 * 0: k > 0 *
46 * *
47 *...........................................................................*
48 * *
49 * distr: Burr distribution Type III *
50 * first parameter is 3 *
51 * *
52 * cdf: F(x) = (x^(-c)+1)^(-k) *
53 * pdf: f(x) = c * k * x^(-c-1) * (x^(-c) + 1)^(-k-1) *
54 * domain: 0 < x < infinity *
55 * constant: 1 *
56 * *
57 * parameters: 2 *
58 * 0: k > 0 *
59 * 1: c > 0 *
60 * *
61 *...........................................................................*
62 * *
63 * distr: Burr distribution Type IV *
64 * first parameter is 4 *
65 * *
66 * cdf: F(x) = (((c-x)/x)^(1/c) + 1)^(-k) *
67 * pdf: f(x) = *
68 * domain: 0 < x < c *
69 * constant: 1 *
70 * *
71 * parameters: 2 *
72 * 0: k > 0 *
73 * 1: c > 0 *
74 * *
75 *...........................................................................*
76 * *
77 * distr: Burr distribution Type V *
78 * first parameter is 5 *
79 * *
80 * cdf: F(x) = (c * exp(-tan(x)) + 1)^(-k) *
81 * pdf: f(x) = *
82 * domain: -pi/2 < x < pi/2 *
83 * constant: *
84 * *
85 * parameters: 2 *
86 * 0: k > 0 *
87 * 1: c > 0 *
88 * *
89 *...........................................................................*
90 * *
91 * distr: Burr distribution Type VI *
92 * first parameter is 6 *
93 * *
94 * cdf: F(x) = (c * exp(-k*sinh(x)) + 1)^(-k) *
95 * pdf: f(x) = *
96 * domain: -infinity < x < infinity *
97 * constant: *
98 * *
99 * parameters: 2 *
100 * 0: k > 0 *
101 * 1: c > 0 *
102 * *
103 *...........................................................................*
104 * *
105 * distr: Burr distribution Type VII *
106 * first parameter is 7 *
107 * *
108 * cdf: F(x) = 2^(-k) * (1 + tanh(x))^k *
109 * pdf: f(x) = *
110 * domain: -infinity < x < infinity *
111 * constant: *
112 * *
113 * parameters: *
114 * 0: k > 0 *
115 * 1: c > 0 *
116 * *
117 *...........................................................................*
118 * *
119 * distr: Burr distribution Type VIII *
120 * first parameter is 8 *
121 * *
122 * cdf: F(x) = (2/pi * arctan(exp(x)))^k *
123 * pdf: f(x) = *
124 * domain: -infinity < x < infinity *
125 * constant: *
126 * *
127 * parameters: 1 *
128 * 0: k > 0 *
129 * *
130 *...........................................................................*
131 * *
132 * distr: Burr distribution Type IX *
133 * first parameter is 9 *
134 * *
135 * cdf: F(x) = 1 - 2 / (2 + c * ((1+exp(x))^k - 1)) *
136 * pdf: f(x) = *
137 * domain: -infinity < x < infinity *
138 * constant: *
139 * *
140 * parameters: 2 *
141 * 0: k > 0 *
142 * 1: c > 0 *
143 * *
144 *...........................................................................*
145 * *
146 * distr: Burr distribution Type X *
147 * first parameter is 10 *
148 * *
149 * cdf: F(x) = (1 - exp(-x^2))^k *
150 * pdf: f(x) = *
151 * domain: 0 < x < infinity *
152 * constant: *
153 * *
154 * parameters: 1 *
155 * 0: k > 0 *
156 * *
157 *...........................................................................*
158 * *
159 * distr: Burr distribution Type XI *
160 * first parameter is 11 *
161 * *
162 * cdf: F(x) = (x - 1/(2*pi) * sin( 2*pi*x))^k *
163 * pdf: f(x) = *
164 * domain: 0 < x < 1 *
165 * constant: *
166 * *
167 * parameters: *
168 * 0: k > 0 *
169 * *
170 *...........................................................................*
171 * *
172 * distr: Burr distribution Type XII *
173 * first parameter is 12 *
174 * *
175 * cdf: F(x) = 1 - (1 + x^c)^(-k) *
176 * pdf: f(x) = *
177 * domain: 0 < x < infinity *
178 * constant: *
179 * *
180 * parameters: *
181 * 0: k > 0 *
182 * 1: c > 0 *
183 * *
184 *****************************************************************************
185 * *
186 * Copyright (c) 2000-2006, 2010 Wolfgang Hoermann and Josef Leydold *
187 * Department of Statistics and Mathematics, WU Wien, Austria *
188 * *
189 * This program is free software; you can redistribute it and/or modify *
190 * it under the terms of the GNU General Public License as published by *
191 * the Free Software Foundation; either version 2 of the License, or *
192 * (at your option) any later version. *
193 * *
194 * This program is distributed in the hope that it will be useful, *
195 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
196 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
197 * GNU General Public License for more details. *
198 * *
199 * You should have received a copy of the GNU General Public License *
200 * along with this program; if not, write to the *
201 * Free Software Foundation, Inc., *
202 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
203 * *
204 *****************************************************************************/
205
206 /*---------------------------------------------------------------------------*/
207
208 #include <unur_source.h>
209 #include <distr/distr_source.h>
210 #include <distr/cont.h>
211 #include <specfunct/unur_specfunct_source.h>
212 #include "unur_distributions.h"
213 #include "unur_distributions_source.h"
214 #include "unur_stddistr.h"
215
216 /*---------------------------------------------------------------------------*/
217 static const char distr_name[] = "burr";
218
219 /* parameters */
220 #define burr_type params[0]
221 #define k params[1]
222 #define c params[2]
223
224 #define DISTR distr->data.cont
225 /* #define NORMCONSTANT (distr->data.cont.norm_constant) */
226
227 /* function prototypes */
228 /* static double _unur_pdf_burr(double x, const UNUR_DISTR *distr); */
229 /* static double _unur_dpdf_burr(double x, const UNUR_DISTR *distr); */
230 static double _unur_cdf_burr(double x, const UNUR_DISTR *distr);
231 static double _unur_invcdf_burr(double u, const UNUR_DISTR *distr);
232
233 /* static int _unur_upd_mode_burr( UNUR_DISTR *distr ); */
234 /* static int _unur_upd_area_burr( UNUR_DISTR *distr ); */
235 static int _unur_set_params_burr( UNUR_DISTR *distr, const double *params, int n_params );
236
237 /*---------------------------------------------------------------------------*/
238
239 double
_unur_cdf_burr(double x,const UNUR_DISTR * distr)240 _unur_cdf_burr( double x, const UNUR_DISTR *distr )
241 {
242 register const double *params = DISTR.params;
243
244 switch ((int) (burr_type + 0.5)) {
245
246 case 1: /* Type I: F(x) = x */
247 if (x<=0.)
248 return 0.;
249 if (x>=1.)
250 return 1.;
251 return x;
252
253 case 2: /* Type II: F(x) = (exp(-x) + 1)^(-k) */
254 return pow(exp(-x) + 1., -k);
255
256 case 3: /* Type III: F(x) = (x^(-c)+1)^(-k) */
257 if (x<=0.)
258 return 0.;
259 return pow( pow(x, -c) + 1., -k);
260
261 case 4: /* Type IV: F(x) = (((c-x)/x)^(1/c) + 1)^(-k) */
262 if (x<=0.)
263 return 0.;
264 if (x>=c)
265 return 1.;
266 return pow( pow( (c-x)/x, 1/c ) + 1., -k);
267
268 case 5: /* Type V: F(x) = (c * exp(-tan(x)) + 1)^(-k) */
269 if (x<=-M_PI/2.)
270 return 0.;
271 if (x>= M_PI/2.)
272 return 1.;
273 return pow( c * exp(-tan(x)) + 1., -k );
274
275 case 6: /* Type VI: F(x) = (c * exp(-k*sinh(x)) + 1)^(-k) */
276 return pow( c * exp(-k*sinh(x)) + 1., -k );
277
278 case 7: /* Type VII: F(x) = 2^(-k) * (1 + tanh(x))^k */
279 return pow( (1. + tanh(x))/2, k );
280
281 case 8: /* Type VIII:F(x) = (2/pi * arctan(exp(x)))^k */
282 return pow( 2./M_PI * atan(exp(x)), k );
283
284 case 9: /* Type IX: F(x) = 1 - 2 / (2 + c * ((1+exp(x))^k - 1)) */
285 return (1. - 2. / (2. + c * (pow(1.+exp(x), k) - 1.)));
286
287 case 10: /* Type X: F(x) = (1 - exp(-x^2))^k */
288 if (x<=0.) return 0.;
289 return pow( 1. - exp(-x*x), k );
290
291 case 11: /* Type XI: F(x) = (x - 1/(2*pi) * sin(2*pi*x))^k */
292 if (x<=0.)
293 return 0.;
294 if (x>=1.)
295 return 1.;
296 return pow( x - 1./(2.*M_PI) * sin( 2. * M_PI * x), k );
297
298 case 12: /* Type XII: F(x) = 1 - (1 + x^c)^(-k) */
299 if (x<=0.)
300 return 0.;
301 return (1. - pow( 1 + pow(x, c), -k ) );
302
303 default:
304 _unur_error(distr_name,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
305 return INFINITY;
306 }
307
308 } /* end of _unur_cdf_burr() */
309
310 /*---------------------------------------------------------------------------*/
311
312 double
_unur_invcdf_burr(double U,const UNUR_DISTR * distr)313 _unur_invcdf_burr( double U, const UNUR_DISTR *distr )
314 {
315 register const double *params = DISTR.params;
316 double Y;
317
318 switch (distr->id) {
319
320 case UNUR_DISTR_BURR_I:
321 return U;
322
323 case UNUR_DISTR_BURR_II:
324 Y = exp( -log(U)/k ); /* U^(-1/k) */
325 return ( -log( Y - 1. ) );
326
327 case UNUR_DISTR_BURR_III:
328 Y = exp( -log(U)/k ); /* U^(-1/k) */
329 return ( exp( -log( Y - 1. )/c ) );
330
331 case UNUR_DISTR_BURR_IV:
332 Y = exp( -log(U)/k ); /* U^(-1/k) */
333 Y = exp( c * log( Y - 1. )) + 1.;
334 return (c/Y);
335
336 case UNUR_DISTR_BURR_V:
337 Y = exp( -log(U)/k ); /* U^(-1/k) */
338 return atan( -log( (Y - 1.) / c ) );
339
340 case UNUR_DISTR_BURR_VI:
341 Y = exp( -log(U)/k ); /* U^(-1/k) */
342 Y = -log( (Y - 1.) / c)/k;
343 return log( Y + sqrt(Y * Y +1.));
344
345 case UNUR_DISTR_BURR_VII:
346 Y = exp( log(U)/k ); /* U^(1/k) */
347 return ( log(2. * Y / (2. - 2.*Y)) / 2. );
348
349 case UNUR_DISTR_BURR_VIII:
350 Y = exp( log(U)/k ); /* U^(1/k) */
351 return ( log( tan( Y * M_PI/2. ) ) );
352
353 case UNUR_DISTR_BURR_IX:
354 Y = 1. + 2. * U / (c * (1.-U));
355 return log( exp( log(Y) / k) - 1. );
356
357 case UNUR_DISTR_BURR_X:
358 Y = exp( log(U)/k ); /* U^(1/k) */
359 return ( sqrt( -log( 1. - Y ) ) );
360
361 case UNUR_DISTR_BURR_XII:
362 Y = exp( -log(1-U)/k ); /* (1-U)^(-1/k) */
363 return ( exp( log( Y - 1.) / c) );
364
365 case UNUR_DISTR_BURR_XI:
366 default:
367 _unur_error(distr_name,UNUR_ERR_SHOULD_NOT_HAPPEN,"");
368 return INFINITY;
369 }
370
371 } /* end of _unur_invcdf_burr() */
372
373 /*---------------------------------------------------------------------------*/
374
375 int
_unur_set_params_burr(UNUR_DISTR * distr,const double * params,int n_params)376 _unur_set_params_burr( UNUR_DISTR *distr, const double *params, int n_params )
377 {
378 CHECK_NULL(params,UNUR_ERR_NULL);
379
380 /* check number of parameters for == 3 */
381 switch (distr->id) {
382 case UNUR_DISTR_BURR_I:
383 if (n_params > 1) {
384 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
385 n_params = 1;
386 }
387 break;
388
389 case UNUR_DISTR_BURR_II:
390 case UNUR_DISTR_BURR_VII:
391 case UNUR_DISTR_BURR_VIII:
392 case UNUR_DISTR_BURR_X:
393 case UNUR_DISTR_BURR_XI:
394 if (n_params < 2) {
395 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few");
396 return UNUR_ERR_DISTR_NPARAMS;
397 }
398 if (n_params > 2) {
399 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
400 n_params = 2;
401 }
402 break;
403
404 case UNUR_DISTR_BURR_III:
405 case UNUR_DISTR_BURR_IV:
406 case UNUR_DISTR_BURR_V:
407 case UNUR_DISTR_BURR_VI:
408 case UNUR_DISTR_BURR_IX:
409 case UNUR_DISTR_BURR_XII:
410 if (n_params < 3) {
411 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few");
412 return UNUR_ERR_DISTR_NPARAMS;
413 }
414 if (n_params > 3) {
415 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
416 n_params = 3;
417 }
418 break;
419
420 default:
421 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"unkown type");
422 return UNUR_ERR_DISTR_NPARAMS;
423 }
424
425 /* check parameters */
426 if (k <= 0. || (c <= 0. && n_params == 3) ) {
427 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"k <= 0 || c <= 0");
428 return UNUR_ERR_DISTR_DOMAIN;
429 }
430
431 /* copy parameters */
432 DISTR.burr_type = burr_type;
433 switch (n_params) {
434 case 3:
435 DISTR.c = c;
436 case 2:
437 DISTR.k = k;
438 default:
439 break;
440 }
441
442 /* number of arguments */
443 DISTR.n_params = n_params;
444
445 /* set (standard) domain */
446 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
447 DISTR.domain[0] = -INFINITY; /* left boundary */
448 DISTR.domain[1] = INFINITY; /* right boundary */
449
450 switch (distr->id) {
451 case UNUR_DISTR_BURR_I:
452 DISTR.domain[0] = 0.; /* left boundary */
453 DISTR.domain[1] = 1.; /* right boundary */
454 break;
455 case UNUR_DISTR_BURR_III:
456 DISTR.domain[0] = 0.; /* left boundary */
457 break;
458 case UNUR_DISTR_BURR_IV:
459 DISTR.domain[0] = 0.; /* left boundary */
460 DISTR.domain[1] = DISTR.c; /* right boundary */
461 break;
462 case UNUR_DISTR_BURR_V:
463 DISTR.domain[0] = -M_PI/2.; /* left boundary */
464 DISTR.domain[1] = M_PI/2.; /* right boundary */
465 break;
466 case UNUR_DISTR_BURR_X:
467 DISTR.domain[0] = 0.; /* left boundary */
468 break;
469 case UNUR_DISTR_BURR_XI:
470 DISTR.domain[0] = 0.; /* left boundary */
471 DISTR.domain[1] = 1.; /* right boundary */
472 break;
473 case UNUR_DISTR_BURR_XII:
474 DISTR.domain[0] = 0.; /* left boundary */
475 break;
476 }
477 }
478
479 /* pointer to inverse CDF */
480 DISTR.invcdf = ( (distr->id != UNUR_DISTR_BURR_XI)
481 ? _unur_invcdf_burr : NULL );
482
483 return UNUR_SUCCESS;
484 } /* end of _unur_set_params_burr() */
485
486 /*---------------------------------------------------------------------------*/
487
488 struct unur_distr *
unur_distr_burr(const double * params,int n_params)489 unur_distr_burr( const double *params, int n_params )
490 {
491 register struct unur_distr *distr;
492
493 if (n_params < 1) {
494 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few");
495 return NULL;
496 }
497
498 /* get new (empty) distribution object */
499 distr = unur_distr_cont_new();
500
501 /* get type of distribution and
502 set distribution id */
503 switch ((int) (burr_type + 0.5)) {
504 case 1: distr->id = UNUR_DISTR_BURR_I; break;
505 case 2: distr->id = UNUR_DISTR_BURR_II; break;
506 case 3: distr->id = UNUR_DISTR_BURR_III; break;
507 case 4: distr->id = UNUR_DISTR_BURR_IV; break;
508 case 5: distr->id = UNUR_DISTR_BURR_V; break;
509 case 6: distr->id = UNUR_DISTR_BURR_VI; break;
510 case 7: distr->id = UNUR_DISTR_BURR_VII; break;
511 case 8: distr->id = UNUR_DISTR_BURR_VIII; break;
512 case 9: distr->id = UNUR_DISTR_BURR_IX; break;
513 case 10: distr->id = UNUR_DISTR_BURR_X; break;
514 case 11: distr->id = UNUR_DISTR_BURR_XI; break;
515 case 12: distr->id = UNUR_DISTR_BURR_XII; break;
516 default:
517 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"type < 1 || type > 12");
518 free( distr ); return NULL;
519 }
520
521 /* name of distribution */
522 distr->name = distr_name;
523
524 /* how to get special generators */
525 /* DISTR.init = _unur_stdgen_burr_init; */
526
527 /* functions */
528 /* DISTR.pdf = _unur_pdf_burr; pointer to PDF */
529 /* DISTR.dpdf = _unur_dpdf_burr; pointer to derivative of PDF */
530 DISTR.cdf = _unur_cdf_burr; /* pointer to CDF */
531 DISTR.invcdf = _unur_invcdf_burr; /* pointer to inverse CDF */
532
533 /* indicate which parameters are set */
534 distr->set = ( UNUR_DISTR_SET_DOMAIN |
535 UNUR_DISTR_SET_STDDOMAIN );
536 /* UNUR_DISTR_SET_MODE | */
537 /* UNUR_DISTR_SET_PDFAREA ); */
538
539 /* set parameters for distribution */
540 if (_unur_set_params_burr(distr,params,n_params)!=UNUR_SUCCESS) {
541 free(distr);
542 return NULL;
543 }
544
545 /* normalization constant: none */
546
547 /* mode and area below p.d.f. */
548 /* DISTR.mode = 0.; */
549 /* DISTR.area = 1.; */
550
551 /* function for setting parameters and updating domain */
552 DISTR.set_params = _unur_set_params_burr;
553
554 /* function for updating derived parameters */
555 /* DISTR.upd_mode = _unur_upd_mode_burr; funct for computing mode */
556 /* DISTR.upd_area = _unur_upd_area_burr; funct for computing area */
557
558 /* return pointer to object */
559 return distr;
560
561 } /* end of unur_distr_burr() */
562
563 /*---------------------------------------------------------------------------*/
564 #undef burr_type
565 #undef k
566 #undef c
567 #undef DISTR
568 /*---------------------------------------------------------------------------*/
569