diff options
Diffstat (limited to 'gsl-1.9/randist/dirichlet.c')
-rw-r--r-- | gsl-1.9/randist/dirichlet.c | 103 |
1 files changed, 103 insertions, 0 deletions
diff --git a/gsl-1.9/randist/dirichlet.c b/gsl-1.9/randist/dirichlet.c new file mode 100644 index 0000000..ad136e1 --- /dev/null +++ b/gsl-1.9/randist/dirichlet.c @@ -0,0 +1,103 @@ +/* randist/dirichlet.c + * + * This program is free software; you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation; either version 2 of the License, or (at + * your option) any later version. + * + * This program is distributed in the hope that it will be useful, but + * WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with this program; if not, write to the Free Software + * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA. + */ + +#include <config.h> +#include <math.h> +#include <gsl/gsl_math.h> +#include <gsl/gsl_rng.h> +#include <gsl/gsl_randist.h> +#include <gsl/gsl_sf_gamma.h> + + +/* The Dirichlet probability distribution of order K-1 is + + p(\theta_1,...,\theta_K) d\theta_1 ... d\theta_K = + (1/Z) \prod_i=1,K \theta_i^{alpha_i - 1} \delta(1 -\sum_i=1,K \theta_i) + + The normalization factor Z can be expressed in terms of gamma functions: + + Z = {\prod_i=1,K \Gamma(\alpha_i)} / {\Gamma( \sum_i=1,K \alpha_i)} + + The K constants, \alpha_1,...,\alpha_K, must be positive. The K parameters, + \theta_1,...,\theta_K are nonnegative and sum to 1. + + The random variates are generated by sampling K values from gamma + distributions with parameters a=\alpha_i, b=1, and renormalizing. + See A.M. Law, W.D. Kelton, Simulation Modeling and Analysis (1991). + + Gavin E. Crooks <gec@compbio.berkeley.edu> (2002) +*/ + +void +gsl_ran_dirichlet (const gsl_rng * r, const size_t K, + const double alpha[], double theta[]) +{ + size_t i; + double norm = 0.0; + + for (i = 0; i < K; i++) + { + theta[i] = gsl_ran_gamma (r, alpha[i], 1.0); + } + + for (i = 0; i < K; i++) + { + norm += theta[i]; + } + + for (i = 0; i < K; i++) + { + theta[i] /= norm; + } +} + + +double +gsl_ran_dirichlet_pdf (const size_t K, + const double alpha[], const double theta[]) +{ + return exp (gsl_ran_dirichlet_lnpdf (K, alpha, theta)); +} + +double +gsl_ran_dirichlet_lnpdf (const size_t K, + const double alpha[], const double theta[]) +{ + /*We calculate the log of the pdf to minimize the possibility of overflow */ + size_t i; + double log_p = 0.0; + double sum_alpha = 0.0; + + for (i = 0; i < K; i++) + { + log_p += (alpha[i] - 1.0) * log (theta[i]); + } + + for (i = 0; i < K; i++) + { + sum_alpha += alpha[i]; + } + + log_p += gsl_sf_lngamma (sum_alpha); + + for (i = 0; i < K; i++) + { + log_p -= gsl_sf_lngamma (alpha[i]); + } + + return log_p; +} |