Dirchlet RNG

["Gavin Crooks" <gec at threeplusone dot com>]

Hi,

I have written a Dirchlet random number generator in the
style of gsl, which I thought might be of some use.


Gavin Crooks

gec@compbio.berkeley.edu
http://threeplusone.com/

Computational Biology Research Group
University of California, Berkeley

--

/* randist/dirchlet.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., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/*#include <config.h> Is this needed?*/
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_sf.h>

/* The Dirchlet 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.

   This code generates a Dirchlet by sampling K values from gamma
   distributions with parameters \alpha_i, and renormalizing. 
   See A.M. Law, W.D. Kelton, 1991, Simulation Modeling and Analysis

   Gavin E. Crooks <gec@compbio.berkeley.edu> (2002)
*/

void
gsl_ran_dirchlet (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_dirchlet_pdf (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;
  double log_Z  =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 exp(log_p);
}

Attachment: dirchlet.c
Description: Binary data

Attachment: test_dirchlet.c
Description: Binary data