Thread: Initializing a gaussian deviate function with a randomly generated number

  1. #1
    Registered User
    Join Date
    Mar 2005
    Posts
    12

    Question Initializing a gaussian deviate function with a randomly generated number

    Ive already had a lot of help on this (thanks especially to Dave Evans) but im getting myself confused, so i thought id re-state my problem for clarity. Apologies for any repitition and many thanks in advance to anyone that looks at this post!

    Overall: Im trying to use a linear random number generator (called ran2, that gives numbers between 0and1) to initialize a gaussian deviate function that uses the Box-Muller method. I want to generate these "gaussian numbers" with a distribution varying between -1 and +1.

    (Might be worth noting that in the end i want to obtain a set of 8 linear random numbers between 0and1 and 8random numbers following a gaussian distribution that ranges between -1 and +1. The two sets of numbers will ideally be "independently random" of each other, i.e the gaussian is triggered by a "nth" linear random number not included in the 8i will use, but this something for the future!)

    The ran2 code works great and gives my 8 very useful numbers.
    I want to initialize this function using one of the numbers from ran2, but i dont know how to link the two things together.. Any advice on how do this would be great- im really out of my depth here!

    Code:
    #include <math.h>
    
    float gasdev(long *idum)
    {
    	float ran2(long *idum);
    	static int iset=0;
    	static float gset;
    	float fac,rsq,v1,v2;
    
    	if  (iset == 0) {
    		do {
    			v1=2.0*ran2(idum)-1.0;
    			v2=2.0*ran2(idum)-1.0;
    			rsq=v1*v1+v2*v2;
    		} while (rsq >= 1.0 || rsq == 0.0);
    		fac=sqrt(-2.0*log(rsq)/rsq);
    		gset=v1*fac;
    		iset=1;
    		return v2*fac;
    	} else {
    		iset=0;
    		return gset;
    	}
    }
    Thanks for looking.

    For completeness here is my running ran2 code (from Numerical Recipes in C and with much help from the boards here!)

    Code:
    #include <stdio.h>
    #include <stdlib.h>
    #include <math.h>
    #include <time.h>
    
    #define IM1 2147483563
    #define IM2 2147483399
    #define AM (1.0/IM1)
    #define IMM1 (IM1-1)
    #define IA1 40014
    #define IA2 40692
    #define IQ1 53668
    #define IQ2 52774
    #define IR1 12211
    #define IR2 3791
    #define NTAB 32
    #define NDIV (1+IMM1/NTAB)
    #define EPS 1.2e-7
    #define RNMX (1.0-EPS)
    
    /*Long period (>2e10^18) random number generator of L'Ecuyer with Bays-Durham shuffle
    and added safeguards. Returns a random deviate between 0.0 and 1.0 (exlclusive of the endpoint values).
    Call with idum  a negative integer to initialize;
    thereafter do not alter idum between successive deviates in a sequence.
    RNMX should approxiamte the largest floating value that is less than 1. */
    
    float ran2(long *idum);
    
    int main()
    {
    
       float answer;
       long *pointer;        /* don't really need this; just pass &idum each time */
       long idum = -time(0); /* get a value for initialization of ran2()          */
       int i;
       pointer=&idum;
       printf("*pointer = %ld\n", *pointer);
    
       for (i = 1; i < 9; i++) {             /*gives 8numbers; was previously i=0 and <10 */
         answer=ran2(pointer);
         printf ("%2d: %f, *pointer = %ld\n", i, answer, *pointer);    /*prints pointer value used too*/
       }
      printf ("\nPress ENTER to continue.\n");
    
       fflush (stdin);
       (void) getchar();
    
       return (0);
    }
    
    
    /*this is the ran2 function listed in the book, straight from the cd*/
    
    float ran2(long *idum)
    {
    	int j;
    	long k;
    	static long idum2=123456789;
    	static long iy=0;
    	static long iv[NTAB];
    	float temp;
    
    	if (*idum <= 0)
       {
    		if (-(*idum) < 1) *idum=1;
    		else *idum = -(*idum);
    		idum2=(*idum);
    		for (j=NTAB+7;j>=0;j--)
          {
    			k=(*idum)/IQ1;
    			*idum=IA1*(*idum-k*IQ1)-k*IR1;
    			if (*idum < 0) *idum += IM1;
    			if (j < NTAB) iv[j] = *idum;
    		}
    		iy=iv[0];
    	}
    	k=(*idum)/IQ1;
    	*idum=IA1*(*idum-k*IQ1)-k*IR1;
    	if (*idum < 0) *idum += IM1;
    	k=idum2/IQ2;
    	idum2=IA2*(idum2-k*IQ2)-k*IR2;
    	if (idum2 < 0) idum2 += IM2;
    	j=iy/NDIV;
    	iy=iv[j]-idum2;
    	iv[j] = *idum;
    	if (iy < 1) iy += IMM1;
    	if ((temp=AM*iy) > RNMX) return RNMX;
    	else return temp;
    }
    
    #undef IM1
    #undef IM2
    #undef AM
    #undef IMM1
    #undef IA1
    #undef IA2
    #undef IQ1
    #undef IQ2
    #undef IR1
    #undef IR2
    #undef NTAB
    #undef NDIV
    #undef EPS
    #undef RNMX
    Apologies for posting lots of code!

  2. #2
    Registered User
    Join Date
    Mar 2004
    Posts
    536
    Quote Originally Posted by Amoreldor
    Overall: Im trying to use a linear random number generator (called ran2, that gives numbers between 0and1) to initialize a gaussian deviate function that uses the Box-Muller method. I want to generate these "gaussian numbers" with a distribution varying between -1 and +1.
    One more time: the gaussian distribution ranges from - infinity to + infiinity. If you want to get numbers from a normal distribution (mean = 0, variance = 1) you can use gasdev(). gasdev() and other practical gaussian variate generators typically give numbers between -6 and 6 (or somewhere therabouts).

    If you want to get numbers from a normal distribution but reject any that are not between -1 and 1, you could try something like this:

    Code:
    #include <stdio.h>
    #include <stdlib.h>
    #include <time.h>
    
    
    int main()
    {
      float gasdev(long *idum);
      float ran2(long *idum);
    
      int i, j;
    
      float gnumber;
      float gasnum[8]; /* will hold the numbers from gasdev */
      float rnumber;    /* will hold the number from ran2 */
    
      long idum = -time(0); /* or use any negative int */
      int failed = 0;
    
    
    
      printf("Intialization value for gasdev = %ld\n", idum);
      rnumber = ran2(&idum);
    
    
      for (i = 0; i < 8; i++) {
         for (j = 0; j < 10000; j++) {
           gnumber = gasdev(&idum);
           if ((gnumber >= -1 && gnumber <= 1)){
             break;
           }
         }
         if (j == 10000) {
           failed = 1;
           printf("Something is dreadfully wrong: gasdev didn't return\n");
           printf("a number between -1 and 1 in 10000 tries\n");
           break;
         }
         else {
           gasnum[i] = gnumber;
         }
      }
      printf("rnum      = %10.6f\n", rnumber);
    
      if (!failed) {
        for (i = 0; i < 8; i++) {
          printf("gasnum[%d] = %10.6f\n", i, gasnum[i]);
        }
      }
      return (0);
    }
    I called ran2() to get the single uniform variate. Then I used a loop to get the eight gaussian variates. For each time in the loop, I called gasdev() repeatedly until it gave a number within the required range. (I put the error escape clause in, since a program error could cause an infinite loop here, not because I would ever expect 10000 calls to gasdev() to fail to give an acceptable number.)

    Note that you only need to initialize idum once in the program.

    You can compile this along with gasdev() and ran2() that you have already tested and see if it gives you acceptable results. If so, then incorporate it into your program. If not, well you can try to tell us (again) what you really need.


    Regards,

    Dave
    Last edited by Dave Evans; 03-17-2005 at 10:17 AM.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Getting an error with OpenGL: collect2: ld returned 1 exit status
    By Lorgon Jortle in forum C++ Programming
    Replies: 6
    Last Post: 05-08-2009, 08:18 PM
  2. prime number program with function
    By mackieinva in forum C Programming
    Replies: 17
    Last Post: 09-20-2007, 08:36 AM
  3. Replies: 28
    Last Post: 07-16-2006, 11:35 PM
  4. Calling a Thread with a Function Pointer.
    By ScrollMaster in forum Windows Programming
    Replies: 6
    Last Post: 06-10-2006, 08:56 AM
  5. Creating a student grade book-how?
    By Hopelessly confused in forum C Programming
    Replies: 5
    Last Post: 10-03-2002, 08:43 PM