Dirchlet RNG

Dirchlet RNG

[Gavin Crooks]

[-- Attachment #1: Type: text/plain, Size: 2725 bytes --]

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 #2: dirchlet.c --]
[-- Type: application/octet-stream, Size: 2463 bytes --]

/* 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 #3: test_dirchlet.c --]
[-- Type: application/octet-stream, Size: 4811 bytes --]

/* modified randist/test.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 James Theiler, Brian Gough
 * 
 * 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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_rng.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_sf.h>

#define N 100000


void gsl_ran_dirchlet (const gsl_rng * r, const int K, const double alpha[], double theta[]);
double gsl_ran_dirchlet_pdf (const int K, double alpha[], double theta[]);

void test_dirchlet_moments(void);

double test_dirchlet(void);
double test_dirchlet_pdf(double x);
void testPDF (double (*f) (void), double (*pdf)(double), const char *name);


gsl_rng *r_global;

int
main (void)
{
  gsl_ieee_env_setup ();

  gsl_rng_env_setup() ;
  r_global = gsl_rng_alloc (gsl_rng_default);

#define FUNC(x)  test_ ## x,                     "test gsl_ran_" #x
#define FUNC2(x) test_ ## x, test_ ## x ## _pdf, "test gsl_ran_" #x

  testPDF (FUNC2(dirchlet));
  test_dirchlet_moments();

  exit (gsl_test_summary());
}



/* Check that the observed means of the dirchlet variables are
   within reasonable statistical errors of their correct values.
*/

#define DIRCHLET_K 10
void
test_dirchlet_moments(void) 
{
  double alpha[DIRCHLET_K];
  double theta[DIRCHLET_K];
  double theta_sum[DIRCHLET_K];
  
  double alpha_sum =0.0;
  double mean, obs_mean, sd, sigma;
  int status;  
  int k,n;


  for(k=0; k<DIRCHLET_K; k++)
    {
      alpha[k] = gsl_ran_exponential(r_global, 0.1);
      alpha_sum += alpha[k];
      theta_sum[k] = 0.0;
    }

  for(n=0; n<N; n++) 
    {
      gsl_ran_dirchlet(r_global, DIRCHLET_K, alpha, theta);
      for(k=0;k<DIRCHLET_K; k++) 
	theta_sum[k] += theta[k];
    }
  
  for(k=0; k<DIRCHLET_K; k++)
    {
      mean = alpha[k]/ alpha_sum;
      sd   = sqrt( (alpha[k] * (1. - alpha[k]/alpha_sum))/(alpha_sum*(alpha_sum+1.)));
      obs_mean = theta_sum[k] / N;
      sigma = sqrt(N)*fabs(mean-obs_mean)/sd;

      status = (sigma>3.);

      gsl_test (status, "test gsl_ran_dirchlet: mean %g expected, %g observed",
	    mean, obs_mean);
    }
}



double
test_dirchlet (void)
{
  /* This is a bit of a lame test, since when K=2, the dirchlet distribution
     becomes a beta distribution */
  
  int i;
  int K =2;
  double alpha[] = {2.5, 5.0};
  double theta[2];

  gsl_ran_dirchlet ( r_global, K, alpha, theta);

  return theta[0];
}

double
test_dirchlet_pdf (double x)
{
  int i;
  int K =2;
  double alpha[] = {2.5, 5.0};
  double theta[2];

  if( x<=0.0 || x>= 1.0) return 0.0; /* Out of range */

  theta[0]=x;
  theta[1] = 1.0 -x;

  return  gsl_ran_dirchlet_pdf (K, alpha, theta);
}





#define BINS 100

void
testPDF (double (*f) (void), double (*pdf)(double), const char *name)
{
  double count[BINS], p[BINS];
  double a = -5.0, b = +5.0 ;
  double dx = (b - a) / BINS ;
  int i,j,status = 0, status_i =0 ;

  for (i = 0; i < BINS; i++)
    count[i] = 0 ;

  for (i = 0; i < N; i++)
    {
      double r = f ();
      if (r < b && r > a)
	{ 
	  j =  (int)((r - a)/dx) ;
	  count[j]++;
	}
    }
  
  for (i = 0; i < BINS; i++)
    {
      /* Compute an approximation to the integral of p(x) from x to
         x+dx using Simpson's rule */

      double x = a + i * dx ;
#define STEPS 100
      double sum = 0 ;
      
      if (fabs(x) < 1e-10) /* hit the origin exactly */
	x = 0.0 ; 
      
      for (j = 1; j < STEPS; j++)
	sum += pdf(x + j * dx / STEPS) ;

      p[i] =  0.5 * (pdf(x) + 2*sum + pdf(x + dx - 1e-7)) * dx / STEPS ;
    }

  for (i = 0; i < BINS; i++)
    {
      double x = a + i * dx ;
      double d = fabs(count[i] - N*p[i]) ;
      if (p[i] != 0)
	{
	  double s = d / sqrt(N*p[i]) ;
	  status_i = (s > 5) && (d > 1) ;
	}
      else
	{
	  status_i = (count[i] != 0) ;
	}
      status |= status_i ;
      if (status_i) 
	gsl_test (status_i, "%s [%g,%g) (%g/%d=%g observed vs %g expected)", 
		  name, x, x+dx, count[i],N,count[i]/N, p[i]) ;
    }

  if (status == 0)
    gsl_test (status, "%s, sampling against pdf over range [%g,%g) ", 
	      name, a, b) ;
}

Strange output "calling fini: /usr/lib/libgsl.so.0"

[Alexei Meremianin]

Hello,

I used GSL to compute two-dimensional integral
using sequentally the procedure "gsl_integration_qag"

It works but the final output (on stdio) of the program 
is quite intriguing...

28655:
28655:  calling fini: /usr/lib/libgsl.so.0
28655:
28655:
28655:  calling fini: /usr/lib/libgslcblas.so.0
28655:
28655:
28655:  calling fini: /lib/libm.so.6
28655:
28655:
28655:  calling fini: /lib/libc.so.6
28655:

The program was compiled with gcc-2.95.3-52 by the command

$gcc -Wall -O2 -I../include -lm -lgsl -lgslcblas -o h3_integral \
num-integral.c

The system is Intel PIII, SUSE 7.3, gsl ver. 1.0, release 138 according 
to "rpm -qi gsl"

What all that means?

-- 
all the best

Alexei Meremianin

Re: Dirchlet RNG

[Brian Gough]

Gavin Crooks writes:
 > Hi,
 > 
 > I have written a Dirchlet random number generator in the
 > style of gsl, which I thought might be of some use.
 > 

Looks good.  Can you send it as a patch against the current CVS (see
the HACKINg file for details), with the test integrated into the
existing test.c.  An entry for the manual doc/randist.texi would be
useful. Thanks.

Re: Strange output "calling fini: /usr/lib/libgsl.so.0"

[Brian Gough]

Alexei Meremianin writes:
 > Hello,
 >  I used GSL to compute two-dimensional integral using sequentally
 > the procedure "gsl_integration_qag"
 >  It works but the final output (on stdio) of the program is quite
 > intriguing...
 >  28655: 28655: calling fini: /usr/lib/libgsl.so.0
 >  The program was compiled with gcc-2.95.3-52 by the command
 >  $gcc -Wall -O2 -I../include -lm -lgsl -lgslcblas -o h3_integral \
 > num-integral.c
 >  The system is Intel PIII, SUSE 7.3, gsl ver. 1.0, release 138
 > according to "rpm -qi gsl"
 >  What all that means?

I don't know.  It doesn't look specific to GSL though.

regards
Brian Gough

Re: Strange output "calling fini: /usr/lib/libgsl.so.0"

[Alexei Meremianin]

Hello,

The strange output dissapeared when I changed
"*params = bla;" with "params = &blah;"


     > Hello, I used GSL to compute two-dimensional integral using sequentally
     > the procedure "gsl_integration_qag"

     > It works but the final output (on stdio) of the program is quite
     > intriguing...

     > 28655: 28655: calling fini: /usr/lib/libgsl.so.0 28655: 28655: 28655:
     > calling fini: /usr/lib/libgslcblas.so.0 28655: 28655: 28655: calling
     > fini: /lib/libm.so.6 28655: 28655: 28655: calling fini: /lib/libc.so.6
     > 28655:

     > The program was compiled with gcc-2.95.3-52 by the command

     > $gcc -Wall -O2 -I../include -lm -lgsl -lgslcblas -o h3_integral \
     > num-integral.c

     > The system is Intel PIII, SUSE 7.3, gsl ver. 1.0, release 138 according
     > to "rpm -qi gsl"

     > What all that means?

     > -- all the best

     > Alexei Meremianin


-- 
Alexei Meremianin
Theoretical Quantendynamics
Albert-Ludwigs-Universitaet
Hermann-Herder-Str. 3
D-79104 Freiburg, Germany
phone: ++49-761-203-5800 (Office)
       ++49-761-881-2250 (Home)
fax:   ++49-761-203-5883
e-mail: Alexei.Meremianin@physik.uni-freiburg.de