* Dirchlet RNG
@ 2002-08-05 12:10 Gavin Crooks
2002-08-06 8:44 ` Strange output "calling fini: /usr/lib/libgsl.so.0" Alexei Meremianin
2002-08-06 10:43 ` Dirchlet RNG Brian Gough
0 siblings, 2 replies; 5+ messages in thread
From: Gavin Crooks @ 2002-08-05 12:10 UTC (permalink / raw)
To: gsl-discuss
[-- 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) ;
}
^ permalink raw reply [flat|nested] 5+ messages in thread
* Strange output "calling fini: /usr/lib/libgsl.so.0"
2002-08-05 12:10 Dirchlet RNG Gavin Crooks
@ 2002-08-06 8:44 ` Alexei Meremianin
2002-08-08 14:33 ` Brian Gough
2002-08-09 1:42 ` Alexei Meremianin
2002-08-06 10:43 ` Dirchlet RNG Brian Gough
1 sibling, 2 replies; 5+ messages in thread
From: Alexei Meremianin @ 2002-08-06 8:44 UTC (permalink / raw)
To: gsl-discuss
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
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: Dirchlet RNG
2002-08-05 12:10 Dirchlet RNG Gavin Crooks
2002-08-06 8:44 ` Strange output "calling fini: /usr/lib/libgsl.so.0" Alexei Meremianin
@ 2002-08-06 10:43 ` Brian Gough
1 sibling, 0 replies; 5+ messages in thread
From: Brian Gough @ 2002-08-06 10:43 UTC (permalink / raw)
To: Gavin Crooks; +Cc: gsl-discuss
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.
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: Strange output "calling fini: /usr/lib/libgsl.so.0"
2002-08-06 8:44 ` Strange output "calling fini: /usr/lib/libgsl.so.0" Alexei Meremianin
@ 2002-08-08 14:33 ` Brian Gough
2002-08-09 1:42 ` Alexei Meremianin
1 sibling, 0 replies; 5+ messages in thread
From: Brian Gough @ 2002-08-08 14:33 UTC (permalink / raw)
To: Alexei Meremianin; +Cc: gsl-discuss
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
^ permalink raw reply [flat|nested] 5+ messages in thread
* Re: Strange output "calling fini: /usr/lib/libgsl.so.0"
2002-08-06 8:44 ` Strange output "calling fini: /usr/lib/libgsl.so.0" Alexei Meremianin
2002-08-08 14:33 ` Brian Gough
@ 2002-08-09 1:42 ` Alexei Meremianin
1 sibling, 0 replies; 5+ messages in thread
From: Alexei Meremianin @ 2002-08-09 1:42 UTC (permalink / raw)
To: gsl-discuss
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
^ permalink raw reply [flat|nested] 5+ messages in thread
end of thread, other threads:[~2002-08-09 8:42 UTC | newest]
Thread overview: 5+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2002-08-05 12:10 Dirchlet RNG Gavin Crooks
2002-08-06 8:44 ` Strange output "calling fini: /usr/lib/libgsl.so.0" Alexei Meremianin
2002-08-08 14:33 ` Brian Gough
2002-08-09 1:42 ` Alexei Meremianin
2002-08-06 10:43 ` Dirchlet RNG Brian Gough
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).