1 /* randist/beta.c
2 *
3 * Copyright (C) 1996, 1997, 1998, 1999, 2000, 2007 James Theiler, Brian Gough
4 *
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 3 of the License, or (at
8 * your option) any later version.
9 *
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
14 *
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
18 */
19
20 #include <config.h>
21 #include <math.h>
22 #include <gsl/gsl_math.h>
23 #include <gsl/gsl_rng.h>
24 #include <gsl/gsl_randist.h>
25 #include <gsl/gsl_sf_gamma.h>
26
27 /* The beta distribution has the form
28
29 p(x) dx = (Gamma(a + b)/(Gamma(a) Gamma(b))) x^(a-1) (1-x)^(b-1) dx
30
31 The method used here is the one described in Knuth */
32
33 double
gsl_ran_beta(const gsl_rng * r,const double a,const double b)34 gsl_ran_beta (const gsl_rng * r, const double a, const double b)
35 {
36 if ( (a <= 1.0) && (b <= 1.0) )
37 {
38 double U, V, X, Y;
39 while (1)
40 {
41 U = gsl_rng_uniform_pos(r);
42 V = gsl_rng_uniform_pos(r);
43 X = pow(U, 1.0/a);
44 Y = pow(V, 1.0/b);
45 if ((X + Y ) <= 1.0)
46 {
47 if (X + Y > 0)
48 {
49 return X/ (X + Y);
50 }
51 else
52 {
53 double logX = log(U)/a;
54 double logY = log(V)/b;
55 double logM = logX > logY ? logX: logY;
56 logX -= logM;
57 logY -= logM;
58 return exp(logX - log(exp(logX) + exp(logY)));
59 }
60 }
61 }
62 }
63 else
64 {
65 double x1 = gsl_ran_gamma (r, a, 1.0);
66 double x2 = gsl_ran_gamma (r, b, 1.0);
67 return x1 / (x1 + x2);
68 }
69 }
70
71 double
gsl_ran_beta_pdf(const double x,const double a,const double b)72 gsl_ran_beta_pdf (const double x, const double a, const double b)
73 {
74 if (x < 0 || x > 1)
75 {
76 return 0 ;
77 }
78 else
79 {
80 double p;
81
82 double gab = gsl_sf_lngamma (a + b);
83 double ga = gsl_sf_lngamma (a);
84 double gb = gsl_sf_lngamma (b);
85
86 if (x == 0.0 || x == 1.0)
87 {
88 if (a > 1.0 && b > 1.0)
89 {
90 p = 0.0;
91 }
92 else
93 {
94 p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
95 }
96 }
97 else
98 {
99 p = exp (gab - ga - gb + log(x) * (a - 1) + log1p(-x) * (b - 1));
100 }
101
102 return p;
103 }
104 }
105