3 * Copyright (C) 2007 Brian Gough
4 * Copyright (C) 2002 Gavin E. Crooks <gec@compbio.berkeley.edu>
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 3 of the License, or (at
9 * your option) any later version.
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
23 #include <gsl/gsl_math.h>
24 #include <gsl/gsl_rng.h>
25 #include <gsl/gsl_randist.h>
26 #include <gsl/gsl_sf_gamma.h>
29 /* The Dirichlet probability distribution of order K-1 is
31 p(\theta_1,...,\theta_K) d\theta_1 ... d\theta_K =
32 (1/Z) \prod_i=1,K \theta_i^{alpha_i - 1} \delta(1 -\sum_i=1,K \theta_i)
34 The normalization factor Z can be expressed in terms of gamma functions:
36 Z = {\prod_i=1,K \Gamma(\alpha_i)} / {\Gamma( \sum_i=1,K \alpha_i)}
38 The K constants, \alpha_1,...,\alpha_K, must be positive. The K parameters,
39 \theta_1,...,\theta_K are nonnegative and sum to 1.
41 The random variates are generated by sampling K values from gamma
42 distributions with parameters a=\alpha_i, b=1, and renormalizing.
43 See A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991).
45 Gavin E. Crooks <gec@compbio.berkeley.edu> (2002)
48 static void ran_dirichlet_small (const gsl_rng * r, const size_t K, const double alpha[], double theta[]);
51 gsl_ran_dirichlet (const gsl_rng * r, const size_t K,
52 const double alpha[], double theta[])
57 for (i = 0; i < K; i++)
59 theta[i] = gsl_ran_gamma (r, alpha[i], 1.0);
62 for (i = 0; i < K; i++)
67 if (norm < GSL_SQRT_DBL_MIN) /* Handle underflow */
69 ran_dirichlet_small (r, K, alpha, theta);
73 for (i = 0; i < K; i++)
80 /* When the values of alpha[] are small, scale the variates to avoid
81 underflow so that the result is not 0/0. Note that the Dirichlet
82 distribution is defined by a ratio of gamma functions so we can
83 take out an arbitrary factor to keep the values in the range of
87 ran_dirichlet_small (const gsl_rng * r, const size_t K,
88 const double alpha[], double theta[])
91 double norm = 0.0, umax = 0;
93 for (i = 0; i < K; i++)
95 double u = log(gsl_rng_uniform_pos (r)) / alpha[i];
99 if (u > umax || i == 0) {
104 for (i = 0; i < K; i++)
106 theta[i] = exp(theta[i] - umax);
109 for (i = 0; i < K; i++)
111 theta[i] = theta[i] * gsl_ran_gamma (r, alpha[i] + 1.0, 1.0);
114 for (i = 0; i < K; i++)
119 for (i = 0; i < K; i++)
130 gsl_ran_dirichlet_pdf (const size_t K,
131 const double alpha[], const double theta[])
133 return exp (gsl_ran_dirichlet_lnpdf (K, alpha, theta));
137 gsl_ran_dirichlet_lnpdf (const size_t K,
138 const double alpha[], const double theta[])
140 /*We calculate the log of the pdf to minimize the possibility of overflow */
143 double sum_alpha = 0.0;
145 for (i = 0; i < K; i++)
147 log_p += (alpha[i] - 1.0) * log (theta[i]);
150 for (i = 0; i < K; i++)
152 sum_alpha += alpha[i];
155 log_p += gsl_sf_lngamma (sum_alpha);
157 for (i = 0; i < K; i++)
159 log_p -= gsl_sf_lngamma (alpha[i]);