1 /*****************************************************************************
2 * *
3 * UNURAN -- Universal Non-Uniform Random number generator *
4 * *
5 *****************************************************************************
6 * *
7 * FILE: c_beta.c *
8 * *
9 * REFERENCES: *
10 * *
11 * [3] N.L. Johnson, S. Kotz and N. Balakrishnan *
12 * Continuous Univariate Distributions, *
13 * Volume 2, 2nd edition *
14 * John Wiley & Sons, Inc., New York, 1995 *
15 * *
16 *****************************************************************************
17 * *
18 * distr: Beta distribution [3; ch.25, p.210] *
19 * *
20 * pdf: f(x) = (x-a)^(p-1) * (b-x)^(q-1) *
21 * domain: a < x < b *
22 * constant: 1 / ( Beta(p,q) * (b-a)^(p+q-1) ) *
23 * *
24 * parameters: 4 *
25 * 0: p > 0 ... shape *
26 * 1: q > 0 ... shape *
27 * 2: a (0) ... location *
28 * 3: b >a (1) ... location *
29 * *
30 *****************************************************************************
31 * *
32 * standard form *
33 * *
34 * pdf: f(x) = x^(p-1) * (1-x)^(q-1) *
35 * domain: 0 < x < 1 *
36 * constant: 1 / Beta(p,q) *
37 * *
38 * parameters: 2 *
39 * 0: p > 0 ... shape *
40 * 1: q > 0 ... shape *
41 * *
42 * 2: a = 0 *
43 * 3: b = 1 *
44 * *
45 *****************************************************************************
46 * *
47 * Copyright (c) 2000-2010 Wolfgang Hoermann and Josef Leydold *
48 * Department of Statistics and Mathematics, WU Wien, Austria *
49 * *
50 * This program is free software; you can redistribute it and/or modify *
51 * it under the terms of the GNU General Public License as published by *
52 * the Free Software Foundation; either version 2 of the License, or *
53 * (at your option) any later version. *
54 * *
55 * This program is distributed in the hope that it will be useful, *
56 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
57 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
58 * GNU General Public License for more details. *
59 * *
60 * You should have received a copy of the GNU General Public License *
61 * along with this program; if not, write to the *
62 * Free Software Foundation, Inc., *
63 * 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA *
64 * *
65 *****************************************************************************/
66
67 /*---------------------------------------------------------------------------*/
68
69 #include <unur_source.h>
70 #include <distr/distr_source.h>
71 #include <distr/cont.h>
72 #include <specfunct/unur_specfunct_source.h>
73 #include "unur_distributions.h"
74 #include "unur_distributions_source.h"
75 #include "unur_stddistr.h"
76
77 /*---------------------------------------------------------------------------*/
78
79 static const char distr_name[] = "beta";
80
81 /*---------------------------------------------------------------------------*/
82 /* parameters */
83 #define p params[0]
84 #define q params[1]
85 #define a params[2]
86 #define b params[3]
87
88 /*---------------------------------------------------------------------------*/
89
90 #define DISTR distr->data.cont
91 #define LOGNORMCONSTANT (distr->data.cont.norm_constant)
92
93 /*---------------------------------------------------------------------------*/
94
95 /* function prototypes */
96 static double _unur_pdf_beta( double x, const UNUR_DISTR *distr );
97 static double _unur_logpdf_beta( double x, const UNUR_DISTR *distr );
98 static double _unur_dpdf_beta( double x, const UNUR_DISTR *distr );
99 static double _unur_dlogpdf_beta( double x, const UNUR_DISTR *distr );
100 static double _unur_cdf_beta( double x, const UNUR_DISTR *distr );
101 #ifdef _unur_SF_invcdf_beta
102 static double _unur_invcdf_beta( double x, const UNUR_DISTR *distr );
103 #endif
104
105 static int _unur_upd_mode_beta( UNUR_DISTR *distr );
106 static int _unur_upd_area_beta( UNUR_DISTR *distr );
107 inline static double _unur_lognormconstant_beta( const double *params, int n_params );
108 static int _unur_set_params_beta( UNUR_DISTR *distr, const double *params, int n_params );
109
110 /*---------------------------------------------------------------------------*/
111
112 double
_unur_pdf_beta(double x,const UNUR_DISTR * distr)113 _unur_pdf_beta(double x, const UNUR_DISTR *distr)
114 {
115 register const double *params = DISTR.params;
116
117 if (DISTR.n_params > 2)
118 /* standardize */
119 x = (x-a) / (b-a);
120
121 /* standard form */
122
123 if (x > 0. && x < 1.)
124 return exp((p-1.)*log(x) + (q-1.)*log(1.-x) - LOGNORMCONSTANT);
125
126 if ((_unur_iszero(x) && _unur_isone(p)) || (_unur_isone(x) && _unur_isone(q)))
127 return exp(-LOGNORMCONSTANT);
128
129 if ((_unur_iszero(x) && p<1.) || (_unur_isone(x) && q<1.))
130 return INFINITY;
131
132 /* out of support */
133 return 0.;
134
135 } /* end of _unur_pdf_beta() */
136
137 /*---------------------------------------------------------------------------*/
138
139 double
_unur_logpdf_beta(double x,const UNUR_DISTR * distr)140 _unur_logpdf_beta(double x, const UNUR_DISTR *distr)
141 {
142 register const double *params = DISTR.params;
143
144 if (DISTR.n_params > 2)
145 /* standardize */
146 x = (x-a) / (b-a);
147
148 /* standard form */
149
150 if (x > 0. && x < 1.)
151 return ((p-1.)*log(x) + (q-1.)*log(1.-x) - LOGNORMCONSTANT);
152
153 if ((_unur_iszero(x) && _unur_isone(p))
154 || (_unur_isone(x) && _unur_isone(q)))
155 return (-LOGNORMCONSTANT);
156
157 if ((_unur_iszero(x) && p<1.) || (_unur_isone(x) && q<1.))
158 return INFINITY;
159
160 /* out of support */
161 return -INFINITY;
162
163 } /* end of _unur_logpdf_beta() */
164
165 /*---------------------------------------------------------------------------*/
166
167 double
_unur_dpdf_beta(double x,const UNUR_DISTR * distr)168 _unur_dpdf_beta(double x, const UNUR_DISTR *distr)
169 {
170 register const double *params = DISTR.params;
171
172 if (DISTR.n_params > 2) {
173 /* standardize */
174 x = (x-a) / (b-a);
175 }
176
177 /* standard form */
178
179 if (x > 0. && x < 1.)
180 return (exp((p-2.)*log(x) + (q-2.)*log(1.-x) - LOGNORMCONSTANT) * ( (p-1.)*(1.-x) - (q-1.)*x ) / (b-a) );
181
182 if (_unur_iszero(x) && _unur_isone(p))
183 return (1.-q)*exp(-LOGNORMCONSTANT)/(b-a);
184
185 if (_unur_iszero(x) && _unur_isfsame(p,2.))
186 return exp(-LOGNORMCONSTANT)/(b-a);
187
188 if (_unur_iszero(x) && p<2.)
189 return (p>1. ? INFINITY : -INFINITY);
190
191 /* if (_unur_iszero(x) && p>2.) */
192 /* return 0.; */
193
194 if (_unur_isone(x) && _unur_isone(q))
195 return (p-1.)*exp(-LOGNORMCONSTANT)/(b-a);
196
197 if (_unur_isone(x) && _unur_isfsame(q,2.))
198 return -exp(-LOGNORMCONSTANT)/(b-a);
199
200 if (_unur_isone(x) && q<2.)
201 return (q>1. ? -INFINITY : INFINITY);
202
203 /* if (_unur_isone(x) && q>2.) */
204 /* return 0.; */
205
206 /* out of support */
207 return 0.;
208
209 } /* end of _unur_dpdf_beta() */
210
211 /*---------------------------------------------------------------------------*/
212
213 double
_unur_dlogpdf_beta(double x,const UNUR_DISTR * distr)214 _unur_dlogpdf_beta(double x, const UNUR_DISTR *distr)
215 {
216 register const double *params = DISTR.params;
217
218 if (DISTR.n_params > 2) {
219 /* standardize */
220 x = (x-a) / (b-a);
221 }
222
223 /* standard form */
224
225 if (x > 0. && x < 1.)
226 return (((p-1.)/x - (q-1.)/(1.-x)) / (b-a));
227
228 if (_unur_iszero(x) && p<1.)
229 return -INFINITY;
230
231 if (_unur_iszero(x) && _unur_isone(p))
232 return (-(q-1.)/((1.-x)*(b-a)));
233
234 if (_unur_iszero(x) && p>1.)
235 return INFINITY;
236
237 if (_unur_isone(x) && q<1.)
238 return INFINITY;
239
240 if (_unur_isone(x) && _unur_isone(q))
241 return ((p-1.)/(b-a));
242
243 if (_unur_isone(x) && q>1.)
244 return -INFINITY;
245
246 /* out of support */
247 return 0.;
248
249 } /* end of _unur_dlogpdf_beta() */
250
251 /*---------------------------------------------------------------------------*/
252
253 double
_unur_cdf_beta(double x,const UNUR_DISTR * distr)254 _unur_cdf_beta(double x, const UNUR_DISTR *distr)
255 {
256 register const double *params = DISTR.params;
257
258 if (DISTR.n_params > 2)
259 /* standardize */
260 x = (x-a) / (b-a);
261
262 /* standard form */
263
264 /* out of support of p.d.f.? */
265 if (x <= 0.) return 0.;
266 if (x >= 1.) return 1.;
267
268 return _unur_SF_incomplete_beta(x,p,q);
269
270 } /* end of _unur_cdf_beta() */
271
272 /*---------------------------------------------------------------------------*/
273
274 #ifdef _unur_SF_invcdf_beta
275
276 double
_unur_invcdf_beta(double x,const UNUR_DISTR * distr)277 _unur_invcdf_beta(double x, const UNUR_DISTR *distr)
278 {
279 const double *params = DISTR.params;
280
281 if (DISTR.n_params == 2)
282 return _unur_SF_invcdf_beta(x,p,q);
283 else
284 return (a + _unur_SF_invcdf_beta(x,p,q))*(b-a);
285
286 } /* end of _unur_invcdf_beta() */
287 #endif
288
289 /*---------------------------------------------------------------------------*/
290
291 int
_unur_upd_mode_beta(UNUR_DISTR * distr)292 _unur_upd_mode_beta( UNUR_DISTR *distr )
293 {
294 register double *params = DISTR.params;
295
296 if (p <= 1. && q > 1.)
297 DISTR.mode = 0.; /* left limit of domain */
298
299 else if (p > 1. && q <= 1.)
300 DISTR.mode = 1.; /* right limit of domain */
301
302 else if (p > 1. && q > 1.)
303 DISTR.mode = (p - 1.) / (p + q - 2.);
304
305 else {
306 /* p.d.f. is not unimodal */
307 DISTR.mode = INFINITY;
308 return UNUR_ERR_DISTR_PROP;
309 }
310
311 if (DISTR.n_params > 2)
312 DISTR.mode = DISTR.mode * (b - a) + a;
313
314 /* mode must be in domain */
315 if (DISTR.mode < DISTR.domain[0])
316 DISTR.mode = DISTR.domain[0];
317 else if (DISTR.mode > DISTR.domain[1])
318 DISTR.mode = DISTR.domain[1];
319
320 /* o.k. */
321 return UNUR_SUCCESS;
322 } /* end of _unur_upd_mode_beta() */
323
324 /*---------------------------------------------------------------------------*/
325
326 int
_unur_upd_area_beta(UNUR_DISTR * distr)327 _unur_upd_area_beta( UNUR_DISTR *distr )
328 {
329 /* log of normalization constant */
330 LOGNORMCONSTANT = _unur_lognormconstant_beta(DISTR.params,DISTR.n_params);
331
332 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
333 DISTR.area = 1.;
334 return UNUR_SUCCESS;
335 }
336
337 /* else */
338 DISTR.area = ( _unur_cdf_beta( DISTR.domain[1],distr)
339 - _unur_cdf_beta( DISTR.domain[0],distr) );
340 return UNUR_SUCCESS;
341
342 } /* end of _unur_upd_area_beta() */
343
344 /*---------------------------------------------------------------------------*/
345
346 double
_unur_lognormconstant_beta(const double * params,int n_params)347 _unur_lognormconstant_beta(const double *params, int n_params)
348 {
349 if (n_params > 2)
350 /* non-standard form */
351 /* log( Beta(p,q) * (b-a) ) */
352 return (_unur_SF_ln_gamma(p) + _unur_SF_ln_gamma(q) - _unur_SF_ln_gamma(p+q) + log(b-a) );
353
354 else
355 /* standard form */
356 /* log( Beta(p,q) ) */
357 return (_unur_SF_ln_gamma(p) + _unur_SF_ln_gamma(q) - _unur_SF_ln_gamma(p+q));
358
359 } /* end of _unur_lognormconstant_beta() */
360
361 /*---------------------------------------------------------------------------*/
362
363 int
_unur_set_params_beta(UNUR_DISTR * distr,const double * params,int n_params)364 _unur_set_params_beta( UNUR_DISTR *distr, const double *params, int n_params )
365 {
366
367 /* check number of parameters for distribution */
368 if (n_params < 2) {
369 _unur_error(distr_name,UNUR_ERR_DISTR_NPARAMS,"too few"); return UNUR_ERR_DISTR_NPARAMS; }
370 if (n_params == 3) {
371 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"");
372 n_params = 2; }
373 if (n_params > 4) {
374 _unur_warning(distr_name,UNUR_ERR_DISTR_NPARAMS,"too many");
375 n_params = 4; }
376 CHECK_NULL(params,UNUR_ERR_NULL);
377
378
379 /* check parameters p and q */
380 if (p <= 0. || q <= 0.) {
381 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"p <= 0 or q <= 0");
382 return UNUR_ERR_DISTR_DOMAIN;
383 }
384
385 /* check parameters a and b */
386 if (n_params > 2 && a >= b) {
387 _unur_error(distr_name,UNUR_ERR_DISTR_DOMAIN,"a >= b");
388 return UNUR_ERR_DISTR_DOMAIN;
389 }
390
391 /* copy parameters for standard form */
392 DISTR.p = p;
393 DISTR.q = q;
394
395 /* copy optional parameters */
396 if (n_params > 2) {
397 DISTR.a = a;
398 DISTR.b = b;
399 }
400 else { /* or use defaults */
401 DISTR.a = 0.; /* default for a */
402 DISTR.b = 1.; /* default for b */
403 }
404
405 /* store number of parameters */
406 DISTR.n_params = n_params;
407
408 /* set (standard) domain */
409 if (distr->set & UNUR_DISTR_SET_STDDOMAIN) {
410 DISTR.domain[0] = DISTR.a; /* left boundary */
411 DISTR.domain[1] = DISTR.b; /* right boundary */
412 }
413
414 return UNUR_SUCCESS;
415 } /* end of _unur_set_params_beta() */
416
417 /*---------------------------------------------------------------------------*/
418
419 struct unur_distr *
unur_distr_beta(const double * params,int n_params)420 unur_distr_beta( const double *params, int n_params )
421 {
422 register struct unur_distr *distr;
423
424 /* get new (empty) distribution object */
425 distr = unur_distr_cont_new();
426
427 /* set distribution id */
428 distr->id = UNUR_DISTR_BETA;
429
430 /* name of distribution */
431 distr->name = distr_name;
432
433 /* how to get special generators */
434 DISTR.init = _unur_stdgen_beta_init;
435
436 /* functions */
437 DISTR.pdf = _unur_pdf_beta; /* pointer to PDF */
438 DISTR.logpdf = _unur_logpdf_beta; /* pointer to logPDF */
439 DISTR.dpdf = _unur_dpdf_beta; /* pointer to derivative of PDF */
440 DISTR.dlogpdf = _unur_dlogpdf_beta; /* pointer to derivative of logPDF */
441 DISTR.cdf = _unur_cdf_beta; /* pointer to CDF */
442 #ifdef _unur_SF_invcdf_beta
443 DISTR.invcdf = _unur_invcdf_beta; /* pointer to inverse CDF */
444 #endif
445
446 /* indicate which parameters are set */
447 distr->set = ( UNUR_DISTR_SET_DOMAIN |
448 UNUR_DISTR_SET_STDDOMAIN |
449 UNUR_DISTR_SET_PDFAREA |
450 UNUR_DISTR_SET_MODE );
451
452 /* set parameters for distribution */
453 if (_unur_set_params_beta(distr,params,n_params)!=UNUR_SUCCESS) {
454 free(distr);
455 return NULL;
456 }
457
458 /* log of normalization constant */
459 LOGNORMCONSTANT = _unur_lognormconstant_beta(DISTR.params,DISTR.n_params);
460
461 /* mode and area below p.d.f. */
462 _unur_upd_mode_beta( distr );
463 DISTR.area = 1.;
464
465 /* function for setting parameters and updating domain */
466 DISTR.set_params = _unur_set_params_beta;
467
468 /* function for updating derived parameters */
469 DISTR.upd_mode = _unur_upd_mode_beta; /* funct for computing mode */
470 DISTR.upd_area = _unur_upd_area_beta; /* funct for computing area */
471
472 /* return pointer to object */
473 return distr;
474
475 } /* end of unur_distr_beta() */
476
477 /*---------------------------------------------------------------------------*/
478 #undef p
479 #undef q
480 #undef a
481 #undef b
482 #undef DISTR
483 /*---------------------------------------------------------------------------*/
484