TITLE Random number generator
COMMENT
$Id: holt_rnd.mod,v 1.1.1.2 1993/08/07 15:36:29 holt Exp $
Revision 1.1 1993/08/04 14:36:28 holt
Initial revision
These functions generate random numbers with various distributions.
They all call the function holt_random(), which is a random number generator
which decides which of the standard random number generators to call.
The seed value is the time if not explicitly specified.
ENDCOMMENT
NEURON {
SUFFIX random
}
INDEPENDENT {t FROM 0 TO 1 WITH 1 (ms)}
VERBATIM
double holt_random();
void holt_seed(int, int);
double holt_normrand(double, double);
double holt_exprand();
int holt_poisrand(double);
ENDVERBATIM
:
: Returns uniformly distributed random numbers on [0,1).
:
FUNCTION Uniform() {
Uniform = holt_random() : Return a random number.
}
:
: Returns Poisson-distributed integers from a Poisson distribution with
: the specified mean.
:
FUNCTION Poisson(mean) {
VERBATIM
_lPoisson = holt_poisrand(_lmean);
ENDVERBATIM
}
:
: Returns a random number from an exponential distribution.
:
FUNCTION Exp() {
Exp = holt_exprand() : Return a random number.
}
:
: Returns a random number from a normal distribution.
:
FUNCTION Normal(mean, std_dev) {
VERBATIM
_lNormal = holt_normrand(_lmean, _lstd_dev); /* Must do in VERBATIM block or
* else another argument gets
* added in front. */
ENDVERBATIM
}
:
: Seed the random number generator, and specify which algorithm to use:
: Arguments:
: 1) The seed value (an integer).
: 2) The algorithm. 0 = drand48(), 1 = random().
:
PROCEDURE Seed(seedval, algorithm) {
VERBATIM
holt_seed((int)_lseedval, (int)_lalgorithm);
ENDVERBATIM
}
VERBATIM
/* A collection of routines for generating random numbers.
* I abstracted some of these routines from the standard scopmath library
* from Duke university so they could use a different random number generator.
*/
int seed = 0; /* Seed for the random number generator. */
#define DRAND48 0 /* Which random number generator to use. */
#ifndef hpux /* HP/UX doesn't come with random(). */
#define RANDOM 1
#define N_GENERATORS 2 /* How many different generators are implemented. */
#else
#define N_GENERATORS 1
#endif
#define DEFAULT_GENERATOR DRAND48 /* Which is used by default. */
int random_generator; /* Code for which generator to use. */
#ifdef hpux
extern double drand48();
#endif
#ifndef __alpha
extern long random();
extern double drand48();
#endif
/*
* Function to seed the random number generator:
* Arguments:
* 1) The seed. If this is 0, the time is used for the seed.
* 2) Which random number generator to use.
* 0 = drand48()
* 1 = random()
*/
void
holt_seed(int seedval, int gen_code)
{
if (seedval == 0) /* No seed explicitly given? */
seedval=3491;
// seedval = time(0); /* Seed based on the time. uncomment on unix/linux systems */
seed = seedval;
if (gen_code >= N_GENERATORS) /* Illegal parameter? */
gen_code = DEFAULT_GENERATOR; /* Don't give an error message. */
random_generator = gen_code; /* Remember which generator to use. */
switch (random_generator) /* Seed is different depending on which */
{ /* algorithm we use. */
case DRAND48: srand48((long)seedval); break;
// #ifndef hpux // can uncomment these three lines on linux/unix systems if desired
// case RANDOM: srandom(seedval); break;
// #endif
}
}
/*
* Function to return a uniformly distributed random number on [0,1).
* Seeds itself automatically with the time if no seed has been provided.
*/
double
holt_random()
{
if (seed == 0) /* No seed provided yet? */
holt_seed(0, DEFAULT_GENERATOR); /* Seed the random number generator. */
switch (random_generator) /* Use the proper generator: */
{
case DRAND48: return drand48();
// #ifndef hpux
// case RANDOM: return (double)random() / (double)0x7fffffff;
// #endif
default: abort();
}
}
/*--------------------------------------------------------------*/
/* */
/* NORMRAND */
/* */
/* Selects a random number from the normal distribution with */
/* specified mean and standard deviation */
/* */
/* Returns: Double precision number from the distribution */
/* */
/* Calling sequence: normrand(mean, std_dev) */
/* */
/* Arguments */
/* Input: mean, double mean of the distribution */
/* */
/* std_dev, double standard deviation of the */
/* distribution */
/* */
/* Output: arguments are unchanged */
/* */
/* Functions called: random */
/* */
/*--------------------------------------------------------------*/
double
holt_normrand(double mean, double std_dev)
{
double s, v1, v2;
s = 1.0;
while (s >= 1.0)
{
v1 = 2.0 * holt_random() - 1.0;
v2 = 2.0 * holt_random() - 1.0;
s = (v1 * v1) + (v2 * v2);
}
v2 = v1 * sqrt(-2.0 * log(s) / s);
return (v2 * std_dev + mean);
}
/*
* This function returns the number of events in a poisson process
* with mean xm during one iteration. Taken from numerical recipes
* in C.
*/
#define PI 3.141592654
float gammln(float);
int holt_poisrand(double xm) {
static float sq,alxm,g,oldm=(-1.0);
double em,tt,y;
if (xm < 12.0) {
if (xm != oldm) {
oldm=xm;
g=exp(-xm);
}
em = -1;
tt=1.0;
do {
++em;
/* tt *= ran1(idum); */
tt *= holt_random();
} while (tt > g);
} else {
if (xm != oldm) {
oldm=xm;
sq=sqrt(2.0*xm);
alxm=log(xm);
g=xm*alxm-gammln(xm+1.0);
}
do {
do {
/* y=tan(PI*ran1(idum)); */
y=tan(PI*holt_random());
em=sq*y+xm;
} while (em < 0.0);
em=floor(em);
tt=0.9*(1.0+y*y)*exp(em*alxm-gammln(em+1.0)-g);
} while (holt_random() > tt);
}
return (int)em;
}
#undef PI
/* (C) Copr. 1986-92 Numerical Recipes Software #.3. */
float gammln(float xx) {
double x,tmp,ser;
static double cof[6]={76.18009173,-86.50532033,24.01409822,
-1.231739516,0.120858003e-2,-0.536382e-5};
int j;
x=xx-1.0;
tmp=x+5.5;
tmp -= (x+0.5)*log(tmp);
ser=1.0;
for (j=0;j<=5;j++) {
x += 1.0;
ser += cof[j]/x;
}
return -tmp+log(2.50662827465*ser);
}
/****************************************************************/
/* */
/* Abstract: exprand() */
/* */
/* Random sample from the exponential distribution exp(-x) */
/* */
/* Calling sequence: exprand() */
/* */
/* Arguments: */
/* Input: none */
/* */
/* Output: none */
/* */
/* Functions called: holt_random() */
/* */
/* Returns: Double precision number from the distribution */
/* */
/* Files accessed: none */
/* */
/* */
/****************************************************************/
double
holt_exprand()
{
return (-log(holt_random()));
}
ENDVERBATIM