/* probc.c */ #include #include #include #include #include #include "probc.h" /**********************************************************************/ int i_factorial ( int n ) { /**********************************************************************/ /* Purpose: I_FACTORIAL returns N!. Definition: N! = Product ( 1 <= I <= N ) I Modified: 02 February 1999 Author: John Burkardt Parameters: Input, int N, the argument of the factorial function. 0 <= N. Output, int I_FACTORIAL, the factorial of N. */ int fact; int i; /* Check. */ if ( n < 0 ) { printf ( "\n" ); printf ( "I_FACTORIAL - Fatal error!\n" ); printf ( " N < 0.\n" ); return 0; } fact = 1; for ( i = 2; i <= n; i++ ) { fact = fact * i; } return fact; } /**********************************************************************/ float ivec_mean ( int n, int x[] ) { /**********************************************************************/ /* Purpose: IVEC_MEAN returns the mean of an int vector. Modified: 01 May 1999 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int X[N], the vector whose mean is desired. Output, float IVEC_MEAN, the mean, or average, of the vector entries. */ int i; float mean; mean = 0.0; for ( i = 0; i < n; i++ ) { mean = mean + ( float ) x[i]; } mean = mean / ( float ) n; return mean; } /**********************************************************************/ float ivec_variance ( int n, int x[] ) { /**********************************************************************/ /* Purpose: IVEC_VARIANCE returns the variance of an int vector. Modified: 01 May 1999 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, int X[N], the vector whose variance is desired. Output, float IVEC_VARIANCE, the variance of the vector entries. */ int i; float mean; float variance; mean = ivec_mean ( n, x ); variance = 0.0; for ( i = 0; i < n; i++ ) { variance = variance + ( ( float ) x[i] - mean ) * ( ( float ) x[i] - mean ); } if ( n > 1 ) { variance = variance / ( float ) ( n - 1 ); } else { variance = 0.0; } return variance; } /**********************************************************************/ float normal_01_sample ( int *iseed ) { /**********************************************************************/ /* NORMAL_01_SAMPLE returns random values from the standard 0,1 normal distribution. Definition: The standard 0,1 normal distribution has mean 0 and standard deviation 1. Method: The Box-Muller method is used. Modified: 21 April 1999 Author: John Burkardt Parameters: Input/output, int *ISEED, a seed for the random number generator. Output, float NORMAL_01_SAMPLE, a normally distributed random value. */ static int need_new = 1; static float x1 = 0.0; static float x2 = 0.0; float y; if ( need_new == 1 ) { x1 = uniform_01_sample ( iseed ); if ( x1 < 0.0000001 ) { x1 = 0.0000001; } x2 = uniform_01_sample ( iseed ); y = sqrt ( - 2.0 * log ( x1 ) ) * cos ( 2.0 * PI * x2 ); need_new = 0; } else { y = sqrt ( - 2.0 * log ( x1 ) ) * sin ( 2.0 * PI * x2 ); need_new = 1; } return y; } /**********************************************************************/ float poisson_cdf ( int k, float a ) { /**********************************************************************/ /* Purpose: POISSON_CDF evaluates the Poisson CDF. Definition: CDF(K,A) is the probability that the number of events observed in a unit time period will be no greater than K, given that the expected number of events in a unit time period is A. Modified: 28 January 1999 Author: John Burkardt Parameters: Input, int K, the argument of the CDF. Input, float A, the parameter of the PDF. 0 < A. Output, float POISSON_CDF, the value of the CDF. */ float cdf; int i; float last; float new; /* Check. */ if ( a <= 0.0 ) { printf ( "\n" ); printf ( "POISSON_CDF - Fatal error!\n" ); printf ( " A <= 0.\n" ); return 0; } /* Special cases. */ if ( k < 0 ) { cdf = 0.0; return cdf; } /* General case. */ new = exp ( - a ); cdf = new; for ( i = 1; i <= k; i++ ) { last = new; new = last * a / ( float ) i; cdf = cdf + new; } return cdf; } /**********************************************************************/ float poisson_mean ( float a ) { /**********************************************************************/ /* Purpose: POISSON_MEAN returns the mean of the Poisson PDF. Modified: 01 February 1999 Author: John Burkardt Parameters: Input, float A, the parameter of the PDF. 0 < A. Output, float POISSON_MEAN, the mean of the PDF. */ if ( a <= 0.0 ) { printf ( "\n" ); printf ( "POISSON_MEAN - Fatal error!\n" ); printf ( " A <= 0.\n" ); return 0; } return a; } /**********************************************************************/ float poisson_pdf ( int k, float a ) { /**********************************************************************/ /* Purpose: POISSON_PDF evaluates the Poisson PDF. Discussion: PDF(K,A) is the probability that the number of events observed in a unit time period will be K, given the expected number of events in a unit time. The parameter A is the expected number of events per unit time. The Poisson PDF is a discrete version of the exponential PDF. The time interval between two Poisson events is a random variable with the exponential PDF. Modified: 01 February 1999 Author: John Burkardt Parameters: Input, int K, the argument of the PDF. Input, float A, the parameter of the PDF. 0 < A. Output, float POISSON_PDF, the value of the PDF. */ float pdf; /* Check. */ if ( a <= 0.0 ) { printf ( "\n" ); printf ( "POISSON_PDF - Fatal error!\n" ); printf ( " A <= 0.\n" ); return 0; } pdf = exp ( - a ) * pow ( a, ( float ) k ) / r_factorial ( k ); return pdf; } /**********************************************************************/ int poisson_sample ( float a, int *iseed ) { /**********************************************************************/ /* Purpose: POISSON_SAMPLE samples the Poisson PDF. Modified: 02 February 1999 Author: John Burkardt Parameters: Input, float A, the parameter of the PDF. 0 < A. Input/output, int *ISEED, the random number generator seed. Output, int POISSON_SAMPLE, a sample of the PDF. */ float cdf; int i; int KMAX = 100; float last; float new; float sum; float sumold; /* Check. */ if ( a <= 0.0 ) { printf ( "\n" ); printf ( "POISSON_SAMPLE - Fatal error!\n" ); printf ( " A <= 0.\n" ); return 0; } /* Pick a random value of CDF. */ cdf = uniform_01_sample ( iseed ); /* Now simply start at K = 0, and find the first value for which CDF(K-1) <= CDF <= CDF(K). */ sum = 0.0; for ( i = 0; i <= KMAX; i++ ) { sumold = sum; if ( i == 0 ) { new = exp ( - a ); sum = new; } else { last = new; new = last * a / ( float ) i; sum = sum + new; } if ( sumold <= cdf && cdf <= sum ) { return i; } } printf ( "\n" ); printf ( "POISSON_SAMPLE - Warning!\n" ); printf ( " Exceeded KMAX = %d\n", KMAX ); return KMAX; } /**********************************************************************/ float poisson_variance ( float a ) { /**********************************************************************/ /* Purpose: POISSON_VARIANCE returns the variance of the Poisson PDF. Modified: 01 February 1999 Author: John Burkardt Parameters: Input, float A, the parameter of the PDF. 0 < A. Output, float POISSON_VARIANCE, the variance of the PDF. */ /* Check. */ if ( a <= 0.0 ) { printf ( "\n" ); printf ( "POISSON_VARIANCE - Fatal error!\n" ); printf ( " A <= 0.\n" ); return 0; } return a; } /**********************************************************************/ float r_factorial ( int n ) { /**********************************************************************/ /* Purpose: R_FACTORIAL returns N!. Definition: N! = Product ( 1 <= I <= N ) I Modified: 02 February 1999 Author: John Burkardt Parameters: Input, int N, the argument of the factorial function. 0 <= N. Output, float R_FACTORIAL, the factorial of N. */ float fact; int i; /* Check. */ if ( n < 0 ) { printf ( "\n" ); printf ( "R_FACTORIAL - Fatal error!\n" ); printf ( " N < 0.\n" ); return 0; } fact = 1.0; for ( i = 2; i <= n; i++ ) { fact = fact * ( float ) i; } return fact; } /**********************************************************************/ float rvec_mean ( int n, float x[] ) { /**********************************************************************/ /* Purpose: RVEC_MEAN returns the mean of a float vector. Modified: 01 May 1999 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, float X[N], the vector whose mean is desired. Output, float RVEC_MEAN, the mean, or average, of the vector entries. */ int i; float mean; mean = 0.0; for ( i = 0; i < n; i++ ) { mean = mean + x[i]; } mean = mean / ( float ) n; return mean; } /**********************************************************************/ float rvec_variance ( int n, float x[] ) { /**********************************************************************/ /* Purpose: RVEC_VARIANCE returns the variance of a float vector. Modified: 01 May 1999 Author: John Burkardt Parameters: Input, int N, the number of entries in the vector. Input, float X[N], the vector whose variance is desired. Output, float RVEC_VARIANCE, the variance of the vector entries. */ int i; float mean; float variance; mean = rvec_mean ( n, x ); variance = 0.0; for ( i = 0; i < n; i++ ) { variance = variance + ( x[i] - mean ) * ( x[i] - mean ); } if ( n > 1 ) { variance = variance / ( float ) ( n - 1 ); } else { variance = 0.0; } return variance; } /**********************************************************************/ float uniform_01_sample ( int *iseed ) { /**********************************************************************/ /* Purpose: UNIFORM_01_SAMPLE is a random number generator. Formula: ISEED = ISEED * (7**5) mod (2**31 - 1) UNIFORM_01_SAMPLE = ISEED * / ( 2**31 - 1 ) Parameters: Input/output, int *ISEED, the integer "seed" used to generate the output random number, and updated in preparation for the next one. *ISEED should not be zero. Output, float UNIFORM_01_SAMPLE, a random value between 0 and 1. */ /* IA = 7**5 IB = 2**15 IB16 = 2**16 IP = 2**31 - 1 */ static int ia = 16807; static int ib15 = 32768; static int ib16 = 65536; static int ip = 2147483647; int iprhi; int ixhi; int k; int leftlo; int loxa; float temp; /* Don't let ISEED be 0. */ if ( *iseed == 0 ) { *iseed = ip; } /* Get the 15 high order bits of ISEED. */ ixhi = *iseed / ib16; /* Get the 16 low bits of ISEED and form the low product. */ loxa = ( *iseed - ixhi * ib16 ) * ia; /* Get the 15 high order bits of the low product. */ leftlo = loxa / ib16; /* Form the 31 highest bits of the full product. */ iprhi = ixhi * ia + leftlo; /* Get overflow past the 31st bit of full product. */ k = iprhi / ib15; /* Assemble all the parts and presubtract IP. The parentheses are essential. */ *iseed = ( ( ( loxa - leftlo * ib16 ) - ip ) + ( iprhi - k * ib15 ) * ib16 ) + k; /* Add IP back in if necessary. */ if ( *iseed < 0 ) { *iseed = *iseed + ip; } /* Multiply by 1 / (2**31-1). */ temp = ( ( float ) *iseed ) * 4.656612875e-10; return ( temp ); }