gaussdev.c File Reference

normally distributed (Gaussian) deviate with zero mean and unit variance. More...

#include "libtpcmodel.h"

Go to the source code of this file.

Macros
#define	RAND_MAX   32767

Functions
double	gaussdev ()
unsigned int	drandSeed (short int seed)
	Make and optionally set the seed for rand(), drand, drandRange, and drandGaussian().
void	init_gaussdev ()
double	gaussdev2 ()
double	drand ()
int	rand_range (int nr, double *d, double low, double up, int type)

Variables
long int	GAUSSDEV_SEED

Detailed Description

normally distributed (Gaussian) deviate with zero mean and unit variance.

Author

Vesa Oikonen

Definition in file gaussdev.c.

Macro Definition Documentation

RAND_MAX

#define RAND_MAX   32767

Max number for rand()

Definition at line 13 of file gaussdev.c.

Referenced by bootstrap(), drand(), and random_shuffle().

Function Documentation

drand()

double drand

(

)

Alternative function to rand() which returns a double precision floating point number in the range of [0,1].

Returns

Random value in the range [0,1].

See also

rand_range, gaussdev

Definition at line 142 of file gaussdev.c.

{
  double d, s;
  s=1.0/(1.0+RAND_MAX);
  do {
    d = ( ( s*rand() + rand() )*s + rand() ) * s;
  } while(d>=1.0);
  return d;
}

Referenced by gaussdev(), gaussdev2(), rand_range(), tgoRandomParameters(), and tgoRandomParametersST().

drandSeed()

unsigned int drandSeed

(

short int

seed

)

Make and optionally set the seed for rand(), drand, drandRange, and drandGaussian().

Uses microseconds from the computer clock and process ID to reduce the chance of getting the same seed for simultaneously executing program threads and instances.

Returns

Returns the seed for srand().

Parameters

seed	Also sets seed with srand (1) or not (0)

Definition at line 63 of file gaussdev.c.

  {
  unsigned int li;
#ifdef HAVE_TIMESPEC_GET
  struct timespec ts;
  timespec_get(&ts, TIME_UTC);
  li=((ts.tv_sec % 10000)*523 ^ ts.tv_nsec*10) ^ ((getpid() % 1000)*983);
#elif defined HAVE_CLOCK_GETTIME
  struct timespec ts;
  clock_gettime(CLOCK_REALTIME, &ts);
  li=((ts.tv_sec % 10000)*523 ^ ts.tv_nsec*10) ^ ((getpid() % 1000)*983);
#elif defined HAVE_GETTIMEOFDAY
  struct timeval tv;
  gettimeofday(&tv, 0);
  li=((tv.tv_sec % 10000)*523 ^ tv.tv_usec*13) ^ ((getpid() % 1000)*983);
#else
  li=(unsigned int)time(NULL)+(unsigned int)getpid();
#endif
  li+=(unsigned int)rand();
  if(seed) srand(li);
  //printf("seed := %u\n", li); printf("RAND_MAX := %u\n", RAND_MAX);
  return(li);
}

Referenced by tgo().

gaussdev()

double gaussdev

(

)

Applies the polar form of Box-Müller transform to produce pseudo-random numbers with Gaussian (normal) distribution which has a zero mean and standard deviation of one. Box GEP, Muller ME. A note on the generation of random normal deviates. Annals of Mathematical Statistics, Volume 29, Issue 2, 1958, 610-611. Available from JSTOR https://www.jstor.org/

Returns

Returns the pseudo-random number.

See also

gaussdev2, drand

Definition at line 30 of file gaussdev.c.

{
  static int ready=0, first=1;
  static double dev;
  double fac, rsq, a, b;
  if(first) {first=0; init_gaussdev();}
 
  /* If we don't have deviate already, then we'll have to make one */
  if(!ready) {
    do {
      //a = 2.*(double)rand()/(double)RAND_MAX - 1.0; 
      //b = 2.*(double)rand()/(double)RAND_MAX - 1.0; 
      a = 2.*drand() - 1.0; 
      b = 2.*drand() - 1.0; 
      rsq = a*a + b*b; 
    } while (rsq>=1.0 || rsq==0.0);
   
    fac = sqrt(-2.0*log(rsq)/rsq);
    dev=a*fac; ready=1; 
    return(b*fac); 
  } else { /* dev is ready so return it */
    ready=0;
    return(dev);
  }
}

gaussdev2()

double gaussdev2

(

)

Applies the polar form of Box-Müller transform to produce pseudo-random numbers with Gaussian (normal) distribution which has a zero mean and standard deviation of one. This function does never set seed, like gaussdev() does, therefore set seed for random number generator before first calling this routine, for example with srand(time(NULL));

Box GEP, Muller ME. A note on the generation of random normal deviates. Annals of Mathematical Statistics, Volume 29, Issue 2, 1958, 610-611. Available from JSTOR https://www.jstor.org/

Returns

Returns the pseudo-random number.

See also

gaussdev, drand

Definition at line 112 of file gaussdev.c.

{
  static int ready=0;
  static double dev;
  double fac, rsq, a, b;
 
  /* If we don't have deviate already, then we'll have to make one */
  if(!ready) {
    do {
      a = 2.*drand() - 1.0; 
      b = 2.*drand() - 1.0; 
      rsq = a*a + b*b; 
    } while (rsq>=1.0 || rsq==0.0);
   
    fac = sqrt(-2.0*log(rsq)/rsq);
    dev=a*fac; ready=1; 
    return(b*fac); 
  } else { /* dev is ready so return it */
    ready=0;
    return(dev);
  }
}

init_gaussdev()

void init_gaussdev

(

)

Initiate random number generator for gaussdev()

Definition at line 92 of file gaussdev.c.

{
  if(GAUSSDEV_SEED<1L) GAUSSDEV_SEED=893165470L;
  srand(GAUSSDEV_SEED);
}

Referenced by gaussdev().

rand_range()

int rand_range	(	int	nr,
		double *	d,
		double	low,
		double	up,
		int	type )

Fills the double array with random numbers between specified limits. Set seed for random number generator before calling this routine, for example with srand(time(NULL));

Returns

Returns 0 when successful, otherwise <> 0.

Parameters

nr	Nr of values in double array
d	Pointer to allocated double array
low	Lower limit for random values
up	Upper limit for random values
type	Distribution: 0=even, 1=square-root transformation

Definition at line 159 of file gaussdev.c.

  {
  int i;
  double dif, v, stl, stu;
 
  if(nr<1) return 0;
  if(d==NULL || type<0 || type>1) return 1;
 
  dif=up-low; if(dif<0.0) return 2;
  if(dif==0.0) {
    for(i=0; i<nr; i++) d[i]=low;
    return 0;
  }
 
  if(type==0) {
    for(i=0; i<nr; i++) d[i] = drand()*dif + low;
  } else if(type==1) {
    stl=copysign(sqrt(fabs(low)),low); if(!isnormal(stl)) stl=0.0;
    stu=copysign(sqrt(fabs(up)), up); if(!isnormal(stu)) stu=0.0;
    dif=stu-stl;
    for(i=0; i<nr; i++) {v=drand()*dif+stl; d[i]=copysign(v*v, v);}
  }
 
  return 0;
}

Variable Documentation

GAUSSDEV_SEED

long int GAUSSDEV_SEED

Seed for random number generator

Definition at line 7 of file gaussdev.c.

Referenced by init_gaussdev().