# include # include # include # include # include "truncated_normal.h" /******************************************************************************/ int i4_uniform_ab ( int a, int b, int *seed ) /******************************************************************************/ /* Purpose: i4_uniform_ab() returns a scaled pseudorandom I4 between A and B. Discussion: The pseudorandom number should be uniformly distributed between A and B. Licensing: This code is distributed under the MIT license. Modified: 24 May 2012 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Second Edition, Springer, 1987, ISBN: 0387964673, LC: QA76.9.C65.B73. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, December 1986, pages 362-376. Pierre L'Ecuyer, Random Number Generation, in Handbook of Simulation, edited by Jerry Banks, Wiley, 1998, ISBN: 0471134031, LC: T57.62.H37. Peter Lewis, Allen Goodman, James Miller, A Pseudo-Random Number Generator for the System/360, IBM Systems Journal, Volume 8, Number 2, 1969, pages 136-143. Input: int A, B, the limits of the interval. int *SEED, a seed for the random number generator. Output: int *SEED: an updated seed. int I4_UNIFORM_AB, a number between A and B. */ { int c; const int i4_huge = 2147483647; int k; float r; int value; if ( *seed == 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "I4_UNIFORM_AB - Fatal error!\n" ); fprintf ( stderr, " Input value of SEED = 0.\n" ); exit ( 1 ); } /* Guaranteee A <= B. */ if ( b < a ) { c = a; a = b; b = c; } k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + i4_huge; } r = ( float ) ( *seed ) * 4.656612875E-10; /* Scale R to lie between A-0.5 and B+0.5. */ r = ( 1.0 - r ) * ( ( float ) ( a ) - 0.5 ) + r * ( ( float ) ( b ) + 0.5 ); /* Round R to the nearest integer. */ value = round ( r ); /* Guarantee that A <= VALUE <= B. */ if ( value < a ) { value = a; } if ( b < value ) { value = b; } return value; } /******************************************************************************/ double normal_01_cdf ( double x ) /******************************************************************************/ /* Purpose: normal_01_cdf() evaluates the Normal 01 CDF. Licensing: This code is distributed under the MIT license. Modified: 10 February 1999 Author: John Burkardt Reference: A G Adams, Areas Under the Normal Curve, Algorithm 39, Computer j., Volume 12, pages 197-198, 1969. Input: double X, the argument of the CDF. Output: double CDF, the value of the CDF. */ { double a1 = 0.398942280444; double a2 = 0.399903438504; double a3 = 5.75885480458; double a4 = 29.8213557808; double a5 = 2.62433121679; double a6 = 48.6959930692; double a7 = 5.92885724438; double b0 = 0.398942280385; double b1 = 3.8052E-08; double b2 = 1.00000615302; double b3 = 3.98064794E-04; double b4 = 1.98615381364; double b5 = 0.151679116635; double b6 = 5.29330324926; double b7 = 4.8385912808; double b8 = 15.1508972451; double b9 = 0.742380924027; double b10 = 30.789933034; double b11 = 3.99019417011; double cdf; double q; double y; /* |X| <= 1.28. */ if ( fabs ( x ) <= 1.28 ) { y = 0.5 * x * x; q = 0.5 - fabs ( x ) * ( a1 - a2 * y / ( y + a3 - a4 / ( y + a5 + a6 / ( y + a7 ) ) ) ); /* 1.28 < |X| <= 12.7 */ } else if ( fabs ( x ) <= 12.7 ) { y = 0.5 * x * x; q = exp ( - y ) * b0 / ( fabs ( x ) - b1 + b2 / ( fabs ( x ) + b3 + b4 / ( fabs ( x ) - b5 + b6 / ( fabs ( x ) + b7 - b8 / ( fabs ( x ) + b9 + b10 / ( fabs ( x ) + b11 ) ) ) ) ) ); /* 12.7 < |X| */ } else { q = 0.0; } /* Take account of negative X. */ if ( x < 0.0 ) { cdf = q; } else { cdf = 1.0 - q; } return cdf; } /******************************************************************************/ double normal_01_cdf_inv ( double p ) /******************************************************************************/ /* Purpose: normal_01_cdf_inv() inverts the standard normal CDF. Discussion: The result is accurate to about 1 part in 10^16. Licensing: This code is distributed under the MIT license. Modified: 27 December 2004 Author: Original FORTRAN77 version by Michael Wichura. This version by John Burkardt. Reference: Michael Wichura, The Percentage Points of the Normal Distribution, Algorithm AS 241, Applied Statistics, Volume 37, Number 3, pages 477-484, 1988. Input: double P, the value of the cumulative probability densitity function. 0 < P < 1. If P is outside this range, an "infinite" value is returned. Output: double NORMAL_01_CDF_INV, the normal deviate value with the property that the probability of a standard normal deviate being less than or equal to this value is P. */ { double a[8] = { 3.3871328727963666080, 1.3314166789178437745E+02, 1.9715909503065514427E+03, 1.3731693765509461125E+04, 4.5921953931549871457E+04, 6.7265770927008700853E+04, 3.3430575583588128105E+04, 2.5090809287301226727E+03 }; double b[8] = { 1.0, 4.2313330701600911252E+01, 6.8718700749205790830E+02, 5.3941960214247511077E+03, 2.1213794301586595867E+04, 3.9307895800092710610E+04, 2.8729085735721942674E+04, 5.2264952788528545610E+03 }; double c[8] = { 1.42343711074968357734, 4.63033784615654529590, 5.76949722146069140550, 3.64784832476320460504, 1.27045825245236838258, 2.41780725177450611770E-01, 2.27238449892691845833E-02, 7.74545014278341407640E-04 }; double const1 = 0.180625; double const2 = 1.6; double d[8] = { 1.0, 2.05319162663775882187, 1.67638483018380384940, 6.89767334985100004550E-01, 1.48103976427480074590E-01, 1.51986665636164571966E-02, 5.47593808499534494600E-04, 1.05075007164441684324E-09 }; double e[8] = { 6.65790464350110377720, 5.46378491116411436990, 1.78482653991729133580, 2.96560571828504891230E-01, 2.65321895265761230930E-02, 1.24266094738807843860E-03, 2.71155556874348757815E-05, 2.01033439929228813265E-07 }; double f[8] = { 1.0, 5.99832206555887937690E-01, 1.36929880922735805310E-01, 1.48753612908506148525E-02, 7.86869131145613259100E-04, 1.84631831751005468180E-05, 1.42151175831644588870E-07, 2.04426310338993978564E-15 }; double q; double r; double split1 = 0.425; double split2 = 5.0; double value; if ( p <= 0.0 ) { value = - HUGE_VAL; return value; } if ( 1.0 <= p ) { value = HUGE_VAL; return value; } q = p - 0.5; if ( fabs ( q ) <= split1 ) { r = const1 - q * q; value = q * r8poly_value_horner ( 7, a, r ) / r8poly_value_horner ( 7, b, r ); } else { if ( q < 0.0 ) { r = p; } else { r = 1.0 - p; } if ( r <= 0.0 ) { value = HUGE_VAL; } else { r = sqrt ( - log ( r ) ); if ( r <= split2 ) { r = r - const2; value = r8poly_value_horner ( 7, c, r ) / r8poly_value_horner ( 7, d, r ); } else { r = r - split2; value = r8poly_value_horner ( 7, e, r ) / r8poly_value_horner ( 7, f, r ); } } if ( q < 0.0 ) { value = - value; } } return value; } /******************************************************************************/ void normal_01_cdf_values ( int *n_data, double *x, double *fx ) /******************************************************************************/ /* Purpose: normal_01_cdf_values() returns some values of the Normal 01 CDF. Discussion: In Mathematica, the function can be evaluated by: Needs["Statistics`ContinuousDistributions`"] dist = NormalDistribution [ 0, 1 ] CDF [ dist, x ] Licensing: This code is distributed under the MIT license. Modified: 28 August 2004 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 17 static double fx_vec[N_MAX] = { 0.5000000000000000E+00, 0.5398278372770290E+00, 0.5792597094391030E+00, 0.6179114221889526E+00, 0.6554217416103242E+00, 0.6914624612740131E+00, 0.7257468822499270E+00, 0.7580363477769270E+00, 0.7881446014166033E+00, 0.8159398746532405E+00, 0.8413447460685429E+00, 0.9331927987311419E+00, 0.9772498680518208E+00, 0.9937903346742239E+00, 0.9986501019683699E+00, 0.9997673709209645E+00, 0.9999683287581669E+00 }; static double x_vec[N_MAX] = { 0.0000000000000000E+00, 0.1000000000000000E+00, 0.2000000000000000E+00, 0.3000000000000000E+00, 0.4000000000000000E+00, 0.5000000000000000E+00, 0.6000000000000000E+00, 0.7000000000000000E+00, 0.8000000000000000E+00, 0.9000000000000000E+00, 0.1000000000000000E+01, 0.1500000000000000E+01, 0.2000000000000000E+01, 0.2500000000000000E+01, 0.3000000000000000E+01, 0.3500000000000000E+01, 0.4000000000000000E+01 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *x = 0.0; *fx = 0.0; } else { *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double normal_01_mean ( ) /******************************************************************************/ /* Purpose: normal_01_mean() returns the mean of the Normal 01 PDF. Licensing: This code is distributed under the MIT license. Modified: 18 September 2004 Author: John Burkardt Output: double NORMAL_01_MEAN, the mean of the PDF. */ { double mean; mean = 0.0; return mean; } /******************************************************************************/ double normal_01_moment ( int order ) /******************************************************************************/ /* Purpose: normal_01_moment() evaluates moments of the Normal 01 PDF. Licensing: This code is distributed under the MIT license. Modified: 31 August 2013 Author: John Burkardt Input: int ORDER, the order of the moment. 0 <= ORDER. Output: double NORMAL_01_MOMENT, the value of the moment. */ { double value; if ( ( order % 2 ) == 0 ) { value = r8_factorial2 ( order - 1 ); } else { value = 0.0; } return value; } /******************************************************************************/ double normal_01_pdf ( double x ) /******************************************************************************/ /* Purpose: normal_01_pdf() evaluates the Normal 01 PDF. Discussion: The Normal 01 PDF is also called the "Standard Normal" PDF, or the Normal PDF with 0 mean and standard deviation 1. PDF(X) = exp ( - 0.5 * X^2 ) / sqrt ( 2 * PI ) Licensing: This code is distributed under the MIT license. Modified: 18 September 2004 Author: John Burkardt Input: double X, the argument of the PDF. Output: double NORMAL_01_PDF, the value of the PDF. */ { double pdf; const double r8_pi = 3.14159265358979323; pdf = exp ( -0.5 * x * x ) / sqrt ( 2.0 * r8_pi ); return pdf; } /******************************************************************************/ double normal_01_sample ( int *seed ) /******************************************************************************/ /* Purpose: normal_01_sample() samples the standard normal probability distribution. Discussion: The standard normal probability distribution function (PDF) has mean 0 and standard deviation 1. The Box-Muller method is used, which is efficient, but generates two values at a time. Licensing: This code is distributed under the MIT license. Modified: 18 September 2004 Author: John Burkardt Input: int *SEED, a seed for the random number generator. Output: int *SEED: an updated seed. double NORMAL_01_SAMPLE, a normally distributed random value. */ { double r1; double r2; const double r8_pi = 3.14159265358979323; static int used = -1; double x; static double y = 0.0; if ( used == -1 ) { used = 0; } /* If we've used an even number of values so far, generate two more, return one, and save one. */ if ( ( used % 2 )== 0 ) { for ( ; ; ) { r1 = r8_uniform_01 ( seed ); if ( r1 != 0.0 ) { break; } } r2 = r8_uniform_01 ( seed ); x = sqrt ( -2.0 * log ( r1 ) ) * cos ( 2.0 * r8_pi * r2 ); y = sqrt ( -2.0 * log ( r1 ) ) * sin ( 2.0 * r8_pi * r2 ); } else { x = y; } used = used + 1; return x; } /******************************************************************************/ double normal_01_variance ( ) /******************************************************************************/ /* Purpose: normal_01_variance() returns the variance of the Normal 01 PDF. Licensing: This code is distributed under the MIT license. Modified: 18 September 2004 Author: John Burkardt Output: double NORMAL_01_VARIANCE, the variance of the PDF. */ { double variance; variance = 1.0; return variance; } /******************************************************************************/ double normal_ms_cdf ( double x, double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_cdf() evaluates the Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 19 September 2004 Author: John Burkardt Input: double X, the argument of the CDF. double MU, SIGMA, the parameters of the PDF. 0.0 < SIGMA. Output: double NORMAL_MS_CDF, the value of the CDF. */ { double cdf; double y; y = ( x - mu ) / sigma; cdf = normal_01_cdf ( y ); return cdf; } /******************************************************************************/ double normal_ms_cdf_inv ( double cdf, double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_cdf_inv() inverts the Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 19 September 2004 Author: John Burkardt Reference: Algorithm AS 111, Applied Statistics, Volume 26, pages 118-121, 1977. Input: double CDF, the value of the CDF. 0.0 <= CDF <= 1.0. double MU, SIGMA, the parameters of the PDF. 0.0 < SIGMA. Output: double NORMAL_MS_CDF_INV, the corresponding argument. */ { double x; double x2; if ( cdf < 0.0 || 1.0 < cdf ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "NORMAL_MS_CDF_INV - Fatal error!\n" ); fprintf ( stderr, " CDF < 0 or 1 < CDF.\n" ); exit ( 1 ); } x2 = normal_01_cdf_inv ( cdf ); x = mu + sigma * x2; return x; } /******************************************************************************/ double normal_ms_mean ( double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_mean() returns the mean of the Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 19 September 2004 Author: John Burkardt Input: double MU, SIGMA, the parameters of the PDF. 0.0 < SIGMA. Output: double NORMAL_MS_MEAN, the mean of the PDF. */ { double mean; mean = mu; return mean; } /******************************************************************************/ double normal_ms_moment ( int order, double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_moment() evaluates moments of the Normal PDF. Discussion: The formula was posted by John D Cook. Order Moment ----- ------ 0 1 1 mu 2 mu^2 + sigma^2 3 mu^3 + 3 mu sigma^2 4 mu^4 + 6 mu^2 sigma^2 + 3 sigma^4 5 mu^5 + 10 mu^3 sigma^2 + 15 mu sigma^4 6 mu^6 + 15 mu^4 sigma^2 + 45 mu^2 sigma^4 + 15 sigma^6 7 mu^7 + 21 mu^5 sigma^2 + 105 mu^3 sigma^4 + 105 mu sigma^6 8 mu^8 + 28 mu^6 sigma^2 + 210 mu^4 sigma^4 + 420 mu^2 sigma^6 + 105 sigma^8 Licensing: This code is distributed under the MIT license. Modified: 30 August 2013 Author: John Burkardt Input: int ORDER, the order of the moment. 0 <= ORDER. double MU, the mean of the distribution. double SIGMA, the standard deviation of the distribution. Output: double NORMAL_MS_MOMENT, the value of the central moment. */ { int j; int j_hi; double value; j_hi = ( order / 2 ); value = 0.0; for ( j = 0; j <= j_hi; j++ ) { value = value + r8_choose ( order, 2 * j ) * r8_factorial2 ( 2 * j - 1 ) * pow ( mu, order - 2 * j ) * pow ( sigma, 2 * j ); } return value; } /******************************************************************************/ double normal_ms_moment_central ( int order, double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_moment_central() evaluates central moments of the Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 30 August 2013 Author: John Burkardt Input: int ORDER, the order of the moment. 0 <= ORDER. double MU, the mean of the distribution. double SIGMA, the standard deviation of the distribution. Output: double NORMAL_MS_MOMENT_CENTRAL, the value of the central moment. */ { double value; if ( ( order % 2 ) == 0 ) { value = r8_factorial2 ( order - 1 ) * pow ( sigma, order ); } else { value = 0.0; } return value; } /******************************************************************************/ double normal_ms_moment_central_values ( int order, double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_moment_central_values(): moments 0 through 10 of the Normal PDF. Discussion: The formula was posted by John D Cook. Order Moment ----- ------ 0 1 1 0 2 sigma^2 3 0 4 3 sigma^4 5 0 6 15 sigma^6 7 0 8 105 sigma^8 9 0 10 945 sigma^10 Licensing: This code is distributed under the MIT license. Modified: 31 August 2013 Author: John Burkardt Input: int ORDER, the order of the moment. 0 <= ORDER <= 10. double MU, the mean of the distribution. double SIGMA, the standard deviation of the distribution. Output: double NORMAL_MS_MOMENT_CENTRAL_VALUES, the value of the central moment. */ { double value; if ( order == 0 ) { value = 1.0; } else if ( order == 1 ) { value = 0.0; } else if ( order == 2 ) { value = pow ( sigma, 2 ); } else if ( order == 3 ) { value = 0.0; } else if ( order == 4 ) { value = 3.0 * pow ( sigma, 4 ); } else if ( order == 5 ) { value = 0.0; } else if ( order == 6 ) { value = 15.0 * pow ( sigma, 6 ); } else if ( order == 7 ) { value = 0.0; } else if ( order == 8 ) { value = 105.0 * pow ( sigma, 8 ); } else if ( order == 9 ) { value = 0.0; } else if ( order == 10 ) { value = 945.0 * pow ( sigma, 10 ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "NORMAL_MS_MOMENT_CENTRAL_VALUES - Fatal error!\n" ); fprintf ( stderr, " Only ORDERS 0 through 10 are available.\n" ); exit ( 1 ); } return value; } /******************************************************************************/ double normal_ms_moment_values ( int order, double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_moment_values() evaluates moments 0 through 8 of the Normal PDF. Discussion: The formula was posted by John D Cook. Order Moment ----- ------ 0 1 1 mu 2 mu^2 + sigma^2 3 mu^3 + 3 mu sigma^2 4 mu^4 + 6 mu^2 sigma^2 + 3 sigma^4 5 mu^5 + 10 mu^3 sigma^2 + 15 mu sigma^4 6 mu^6 + 15 mu^4 sigma^2 + 45 mu^2 sigma^4 + 15 sigma^6 7 mu^7 + 21 mu^5 sigma^2 + 105 mu^3 sigma^4 + 105 mu sigma^6 8 mu^8 + 28 mu^6 sigma^2 + 210 mu^4 sigma^4 + 420 mu^2 sigma^6 + 105 sigma^8 Licensing: This code is distributed under the MIT license. Modified: 31 August 2013 Author: John Burkardt Input: int ORDER, the order of the moment. 0 <= ORDER <= 8. double MU, the mean of the distribution. double SIGMA, the standard deviation of the distribution. Output: double NORMAL_MS_MOMENT_VALUES, the value of the central moment. */ { double value; if ( order == 0 ) { value = 1.0; } else if ( order == 1 ) { value = mu; } else if ( order == 2 ) { value = pow ( mu, 2 ) + pow ( sigma, 2 ); } else if ( order == 3 ) { value = pow ( mu, 3 ) + 3.0 * mu * pow ( sigma, 2 ); } else if ( order == 4 ) { value = pow ( mu, 4 ) + 6.0 * pow ( mu, 2 ) * pow ( sigma, 2 ) + 3.0 * pow ( sigma, 4 ); } else if ( order == 5 ) { value = pow ( mu, 5 ) + 10.0 * pow ( mu, 3 ) * pow ( sigma, 2 ) + 15.0 * mu * pow ( sigma, 4 ); } else if ( order == 6 ) { value = pow ( mu, 6 ) + 15.0 * pow ( mu, 4 ) * pow ( sigma, 2 ) + 45.0 * pow ( mu, 2 ) * pow ( sigma, 4 ) + 15.0 * pow ( sigma, 6 ); } else if ( order == 7 ) { value = pow ( mu, 7 ) + 21.0 * pow ( mu, 5 ) * pow ( sigma, 2 ) + 105.0 * pow ( mu, 3 ) * pow ( sigma, 4 ) + 105.0 * mu * pow ( sigma, 6 ); } else if ( order == 8 ) { value = pow ( mu, 8 ) + 28.0 * pow ( mu, 6 ) * pow ( sigma, 2 ) + 210.0 * pow ( mu, 4 ) * pow ( sigma, 4 ) + 420.0 * pow ( mu, 2 ) * pow ( sigma, 6 ) + 105.0 * pow ( sigma, 8 ); } else { fprintf ( stderr, "\n" ); fprintf ( stderr, "NORMAL_MS_MOMENT_VALUES - Fatal error!\n" ); fprintf ( stderr, " Only ORDERS 0 through 8 are available.\n" ); exit ( 1 ); } return value; } /******************************************************************************/ double normal_ms_pdf ( double x, double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_pdf() evaluates the Normal PDF. Discussion: PDF(MU,SIGMA;X) = exp ( - 0.5 * ( ( X - MU ) / SIGMA )^2 ) / ( SIGMA * SQRT ( 2 * PI ) ) The normal PDF is also known as the Gaussian PDF. Licensing: This code is distributed under the MIT license. Modified: 19 September 2004 Author: John Burkardt Input: double X, the argument of the PDF. double MU, SIGMA, the parameters of the PDF. 0.0 < SIGMA. Output: double NORMAL_MS_PDF, the value of the PDF. */ { double pdf; const double r8_pi = 3.14159265358979323; double y; y = ( x - mu ) / sigma; pdf = exp ( - 0.5 * y * y ) / ( sigma * sqrt ( 2.0 * r8_pi ) ); return pdf; } /******************************************************************************/ double normal_ms_sample ( double mu, double sigma, int *seed ) /******************************************************************************/ /* Purpose: normal_ms_sample() samples the Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 20 November 2004 Author: John Burkardt Input: double MU, SIGMA, the parameters of the PDF. 0.0 < SIGMA. int *SEED, a seed for the random number generator. Output: int *SEED: an updated seed. double NORMAL_MS_SAMPLE, a sample of the PDF. */ { double x; x = normal_01_sample ( seed ); x = mu + sigma * x; return x; } /******************************************************************************/ double normal_ms_variance ( double mu, double sigma ) /******************************************************************************/ /* Purpose: normal_ms_variance() returns the variance of the Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 19 September 2004 Author: John Burkardt Input: double MU, SIGMA, the parameters of the PDF. 0.0 < SIGMA. Output: double NORMAL_MS_VARIANCE, the variance of the PDF. */ { double variance; variance = sigma * sigma; return variance; } /******************************************************************************/ double r8_choose ( int n, int k ) /******************************************************************************/ /* Purpose: r8_choose() computes the binomial coefficient C(N,K) as an R8. Discussion: The value is calculated in such a way as to avoid overflow and roundoff. The calculation is done in R8 arithmetic. The formula used is: C(N,K) = N! / ( K! * (N-K)! ) Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Reference: ML Wolfson, HV Wright, Algorithm 160: Combinatorial of M Things Taken N at a Time, Communications of the ACM, Volume 6, Number 4, April 1963, page 161. Input: int N, K, the values of N and K. Output: double R8_CHOOSE, the number of combinations of N things taken K at a time. */ { int i; int mn; int mx; double value; if ( k < n - k ) { mn = k; mx = n - k; } else { mn = n - k; mx = k; } if ( mn < 0 ) { value = 0.0; } else if ( mn == 0 ) { value = 1.0; } else { value = ( double ) ( mx + 1 ); for ( i = 2; i <= mn; i++ ) { value = ( value * ( double ) ( mx + i ) ) / ( double ) i; } } return value; } /******************************************************************************/ double r8_factorial2 ( int n ) /******************************************************************************/ /* Purpose: r8_factorial2() computes the double factorial function. Discussion: FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 ) (N even) = Product ( N * (N-2) * (N-4) * ... * 1 ) (N odd) N Value -- ----- 0 1 1 1 2 2 3 3 4 8 5 15 6 48 7 105 8 384 9 945 10 3840 Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Input: int N, the argument of the double factorial function. If N is less than 1, R8_FACTORIAL2 is returned as 1.0. Output: double R8_FACTORIAL2, the value of Factorial2(N). */ { int n_copy; double value; value = 1.0; if ( n < 1 ) { return value; } n_copy = n; while ( 1 < n_copy ) { value = value * ( double ) n_copy; n_copy = n_copy - 2; } return value; } /******************************************************************************/ void r8_factorial2_values ( int *n_data, int *n, double *f ) /******************************************************************************/ /* Purpose: r8_factorial2_values() returns values of the double factorial function. Formula: FACTORIAL2( N ) = Product ( N * (N-2) * (N-4) * ... * 2 ) (N even) = Product ( N * (N-2) * (N-4) * ... * 1 ) (N odd) In Mathematica, the function can be evaluated by: n!! Example: N N!! 0 1 1 1 2 2 3 3 4 8 5 15 6 48 7 105 8 384 9 945 10 3840 Licensing: This code is distributed under the MIT license. Modified: 07 February 2015 Author: John Burkardt Reference: Milton Abramowitz, Irene Stegun, Handbook of Mathematical Functions, National Bureau of Standards, 1964, ISBN: 0-486-61272-4, LC: QA47.A34. Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Daniel Zwillinger, CRC Standard Mathematical Tables and Formulae, 30th Edition, CRC Press, 1996, page 16. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. int *N, the argument of the function. double *F, the value of the function. */ { # define N_MAX 16 static double f_vec[N_MAX] = { 1.0, 1.0, 2.0, 3.0, 8.0, 15.0, 48.0, 105.0, 384.0, 945.0, 3840.0, 10395.0, 46080.0, 135135.0, 645120.0, 2027025.0 }; static int n_vec[N_MAX] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *n = 0; *f = 0.0; } else { *n = n_vec[*n_data-1]; *f = f_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double r8_log_2 ( double x ) /******************************************************************************/ /* Purpose: r8_log_2() returns the logarithm base 2 of an R8. Discussion: R8_LOG_2 ( X ) = Log ( |X| ) / Log ( 2.0 ) Licensing: This code is distributed under the MIT license. Modified: 12 May 2006 Author: John Burkardt Input: double X, the number whose base 2 logarithm is desired. X should not be 0. Output: double R8_LOG_2, the logarithm base 2 of the absolute value of X. It should be true that |X| = 2^R_LOG_2. */ { double value; if ( x == 0.0 ) { value = - HUGE_VAL; } else { value = log ( fabs ( x ) ) / log ( 2.0 ); } return value; } /******************************************************************************/ double r8_mop ( int i ) /******************************************************************************/ /* Purpose: r8_mop() returns the I-th power of -1 as an R8 value. Discussion: An R8 is an double value. Licensing: This code is distributed under the MIT license. Modified: 01 July 2008 Author: John Burkardt Input: int I, the power of -1. Output: double R8_MOP, the I-th power of -1. */ { double value; if ( ( i % 2 ) == 0 ) { value = + 1.0; } else { value = - 1.0; } return value; } /******************************************************************************/ double r8_uniform_01 ( int *seed ) /******************************************************************************/ /* Purpose: r8_uniform_01() returns a unit pseudorandom R8. Discussion: This routine implements the recursion seed = 16807 * seed mod ( 2^31 - 1 ) unif = seed / ( 2^31 - 1 ) The integer arithmetic never requires more than 32 bits, including a sign bit. Licensing: This code is distributed under the MIT license. Modified: 11 August 2004 Author: John Burkardt Reference: Paul Bratley, Bennett Fox, Linus Schrage, A Guide to Simulation, Springer Verlag, pages 201-202, 1983. Bennett Fox, Algorithm 647: Implementation and Relative Efficiency of Quasirandom Sequence Generators, ACM Transactions on Mathematical Software, Volume 12, Number 4, pages 362-376, 1986. Input: int *SEED, a seed for the random number generator. Output: int *SEED: an updated seed. double R8_UNIFORM_01, a new pseudorandom variate, strictly between 0 and 1. */ { int k; double r; k = *seed / 127773; *seed = 16807 * ( *seed - k * 127773 ) - k * 2836; if ( *seed < 0 ) { *seed = *seed + 2147483647; } /* Although SEED can be represented exactly as a 32 bit integer, it generally cannot be represented exactly as a 32 bit real number! */ r = ( double ) ( *seed ) * 4.656612875E-10; return r; } /******************************************************************************/ void r8poly_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8poly_print() prints out a polynomial. Licensing: This code is distributed under the MIT license. Modified: 26 February 2015 Author: John Burkardt Input: int N, the dimension of A. double A[N+1], the polynomial coefficients. A(0) is the constant term and A(N) is the coefficient of X**N. char *TITLE, a title. */ { int i; double mag; char plus_minus; printf ( "\n" ); printf ( "%s\n", title ); printf ( "\n" ); if ( n <= 0 ) { printf ( " p(x) = 0\n" ); return; } if ( a[n] < 0.0 ) { plus_minus = '-'; } else { plus_minus = ' '; } mag = fabs ( a[n] ); if ( 2 <= n ) { printf ( " p(x) = %c%g * x^%d\n", plus_minus, mag, n ); } else if ( n == 1 ) { printf ( " p(x) = %c%g * x\n", plus_minus, mag ); } else if ( n == 0 ) { printf ( " p(x) = %c%g\n", plus_minus, mag ); } for ( i = n - 1; 0 <= i; i-- ) { if ( a[i] < 0.0 ) { plus_minus = '-'; } else { plus_minus = '+'; } mag = fabs ( a[i] ); if ( mag != 0.0 ) { if ( 2 <= i ) { printf ( " %c%g * x^%d\n", plus_minus, mag, i ); } else if ( i == 1 ) { printf ( " %c%g * x\n", plus_minus, mag ); } else if ( i == 0 ) { printf ( " %c%g\n", plus_minus, mag ); } } } return; } /******************************************************************************/ double r8poly_value_horner ( int m, double c[], double x ) /******************************************************************************/ /* Purpose: r8poly_value_horner() evaluates a polynomial using Horner's method. Discussion: The polynomial p(x) = c0 + c1 * x + c2 * x^2 + ... + cm * x^m is to be evaluated at the value X. Licensing: This code is distributed under the MIT license. Modified: 02 January 2015 Author: John Burkardt Input: int M, the degree of the polynomial. double C[M+1], the coefficients of the polynomial. C[0] is the constant term. double X, the point at which the polynomial is to be evaluated. Output: double R8POLY_VALUE_HORNER, the value of the polynomial at X. */ { int i; double value; value = c[m]; for ( i = m - 1; 0 <= i; i-- ) { value = value * x + c[i]; } return value; } /******************************************************************************/ double *r8vec_linspace_new ( int n, double a, double b ) /******************************************************************************/ /* Purpose: r8vec_linspace_new() creates a vector of linearly spaced values. Discussion: An R8VEC is a vector of R8's. 4 points evenly spaced between 0 and 12 will yield 0, 4, 8, 12. In other words, the interval is divided into N-1 even subintervals, and the endpoints of intervals are used as the points. Licensing: This code is distributed under the MIT license. Modified: 29 March 2011 Author: John Burkardt Input: int N, the number of entries in the vector. double A, B, the first and last entries. Output: double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data. */ { int i; double *x; x = ( double * ) malloc ( n * sizeof ( double ) ); if ( n == 1 ) { x[0] = ( a + b ) / 2.0; } else { for ( i = 0; i < n; i++ ) { x[i] = ( ( double ) ( n - 1 - i ) * a + ( double ) ( i ) * b ) / ( double ) ( n - 1 ); } } return x; } /******************************************************************************/ double r8vec_max ( int n, double *dvec ) /******************************************************************************/ /* Purpose: r8vec_max() returns the value of the maximum element in an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 15 October 1998 Author: John Burkardt Input: int N, the number of entries in the array. double *RVEC, a pointer to the first entry of the array. Output: double R8VEC_MAX, the value of the maximum element. This is set to 0.0 if N <= 0. */ { int i; double rmax; double *r8vec_pointer; if ( n <= 0 ) { return 0.0; } for ( i = 0; i < n; i++ ) { if ( i == 0 ) { rmax = *dvec; r8vec_pointer = dvec; } else { r8vec_pointer++; if ( rmax < *r8vec_pointer ) { rmax = *r8vec_pointer; } } } return rmax; } /******************************************************************************/ double r8vec_mean ( int n, double x[] ) /******************************************************************************/ /* Purpose: r8vec_mean() returns the mean of an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 01 May 1999 Author: John Burkardt Input: int N, the number of entries in the vector. double X[N], the vector whose mean is desired. Output: double R8VEC_MEAN, the mean, or average, of the vector entries. */ { int i; double mean; mean = 0.0; for ( i = 0; i < n; i++ ) { mean = mean + x[i]; } mean = mean / ( double ) n; return mean; } /******************************************************************************/ double r8vec_min ( int n, double *dvec ) /******************************************************************************/ /* Purpose: r8vec_min() returns the value of the minimum element in an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 15 October 1998 Author: John Burkardt Input: int N, the number of entries in the array. double *RVEC, a pointer to the first entry of the array. Output: double R8VEC_MIN, the value of the minimum element. This is set to 0.0 if N <= 0. */ { int i; double rmin; double *r8vec_pointer; if ( n <= 0 ) { return 0.0; } for ( i = 0; i < n; i++ ) { if ( i == 0 ) { rmin = *dvec; r8vec_pointer = dvec; } else { r8vec_pointer++; if ( *r8vec_pointer < rmin ) { rmin = *r8vec_pointer; } } } return rmin; } /******************************************************************************/ void r8vec_print ( int n, double a[], char *title ) /******************************************************************************/ /* Purpose: r8vec_print() prints an R8VEC. Discussion: An R8VEC is a vector of R8's. Licensing: This code is distributed under the MIT license. Modified: 08 April 2009 Author: John Burkardt Input: int N, the number of components of the vector. double A[N], the vector to be printed. char *TITLE, a title. */ { int i; fprintf ( stdout, "\n" ); fprintf ( stdout, "%s\n", title ); fprintf ( stdout, "\n" ); for ( i = 0; i < n; i++ ) { fprintf ( stdout, " %8d: %14g\n", i, a[i] ); } return; } /******************************************************************************/ double r8vec_variance ( int n, double x[] ) /******************************************************************************/ /* Purpose: r8vec_variance() returns the variance of an R8VEC. Licensing: This code is distributed under the MIT license. Modified: 01 May 1999 Author: John Burkardt Input: int N, the number of entries in the vector. double X[N], the vector whose variance is desired. Output: double R8VEC_VARIANCE, the variance of the vector entries. */ { int i; double mean; double variance; mean = r8vec_mean ( n, x ); variance = 0.0; for ( i = 0; i < n; i++ ) { variance = variance + ( x[i] - mean ) * ( x[i] - mean ); } if ( 1 < n ) { variance = variance / ( double ) ( n - 1 ); } else { variance = 0.0; } return variance; } /******************************************************************************/ double truncated_normal_ab_cdf ( double x, double mu, double sigma, double a, double b ) /******************************************************************************/ /* Purpose: truncated_normal_ab_cdf() evaluates the truncated Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 24 January 2017 Author: John Burkardt Input: double X, the argument of the CDF. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, B, the lower and upper truncation limits. Output: double TRUNCATED_NORMAL_AB_CDF, the value of the CDF. */ { double alpha; double alpha_cdf; double beta; double beta_cdf; double cdf; double xi; double xi_cdf; if ( x < a ) { cdf = 0.0; } else if ( x <= b ) { alpha = ( a - mu ) / sigma; beta = ( b - mu ) / sigma; xi = ( x - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); beta_cdf = normal_01_cdf ( beta ); xi_cdf = normal_01_cdf ( xi ); cdf = ( xi_cdf - alpha_cdf ) / ( beta_cdf - alpha_cdf ); } else { cdf = 1.0; } return cdf; } /******************************************************************************/ void truncated_normal_ab_cdf_values ( int *n_data, double *mu, double *sigma, double *a, double *b, double *x, double *fx ) /******************************************************************************/ /* Purpose: truncated_normal_ab_cdf_values(): values of the Truncated Normal CDF. Discussion: The Normal distribution, with mean Mu and standard deviation Sigma, is truncated to the interval [A,B]. Licensing: This code is distributed under the MIT license. Modified: 13 September 2013 Author: John Burkardt Reference: Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *MU, the mean of the distribution. double *SIGMA, the standard deviation of the distribution. double *A, *B, the lower and upper truncation limits. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 11 static double a_vec[N_MAX] = { 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0 }; static double b_vec[N_MAX] = { 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0 }; static double fx_vec[N_MAX] = { 0.3371694242213513, 0.3685009225506048, 0.4006444233448185, 0.4334107066903040, 0.4665988676496338, 0.5000000000000000, 0.5334011323503662, 0.5665892933096960, 0.5993555766551815, 0.6314990774493952, 0.6628305757786487 }; static double mu_vec[N_MAX] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 }; static double sigma_vec[N_MAX] = { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 }; static double x_vec[N_MAX] = { 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0, 110.0 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *b = 0.0; *mu = 0.0; *sigma = 0.0; *x = 0.0; *fx = 0.0; } else { *a = a_vec[*n_data-1]; *b = b_vec[*n_data-1]; *mu = mu_vec[*n_data-1]; *sigma = sigma_vec[*n_data-1]; *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double truncated_normal_ab_cdf_inv ( double cdf, double mu, double sigma, double a, double b ) /******************************************************************************/ /* Purpose: truncated_normal_ab_cdf_inv() inverts the truncated Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 14 August 2013 Author: John Burkardt Input: double CDF, the value of the CDF. 0.0D+00 <= CDF <= 1.0. double MU, S, the mean and standard deviation of the parent Normal distribution. double A, B, the lower and upper truncation limits. Output: double TRUNCATED_NORMAL_AB_CDF_INV, the corresponding argument. */ { double alpha; double alpha_cdf; double beta; double beta_cdf; double x; double xi; double xi_cdf; if ( cdf < 0.0 || 1.0 < cdf ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_AB_CDF_INV - Fatal error!\n" ); fprintf ( stderr, " CDF < 0 or 1 < CDF.\n" ); exit ( 1 ); } alpha = ( a - mu ) / sigma; beta = ( b - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); beta_cdf = normal_01_cdf ( beta ); xi_cdf = ( beta_cdf - alpha_cdf ) * cdf + alpha_cdf; xi = normal_01_cdf_inv ( xi_cdf ); x = mu + sigma * xi; return x; } /******************************************************************************/ double truncated_normal_ab_mean ( double mu, double sigma, double a, double b ) /******************************************************************************/ /* Purpose: truncated_normal_ab_mean() returns the mean of the truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 14 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, B, the lower and upper truncation limits. Output: double TRUNCATED_NORMAL_AB_MEAN, the mean of the PDF. */ { double alpha; double alpha_cdf; double alpha_pdf; double beta; double beta_cdf; double beta_pdf; double mean; alpha = ( a - mu ) / sigma; beta = ( b - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); beta_cdf = normal_01_cdf ( beta ); alpha_pdf = normal_01_pdf ( alpha ); beta_pdf = normal_01_pdf ( beta ); mean = mu + sigma * ( alpha_pdf - beta_pdf ) / ( beta_cdf - alpha_cdf ); return mean; } /******************************************************************************/ double truncated_normal_ab_moment ( int order, double mu, double sigma, double a, double b ) /******************************************************************************/ /* Purpose: truncated_normal_ab_moment(): moments of the truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 11 September 2013 Author: John Burkardt Reference: Phoebus Dhrymes, Moments of Truncated Normal Distributions, May 2005. Input: int ORDER, the order of the moment. 0 <= ORDER. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. 0.0 < S. double A, B, the lower and upper truncation limits. A < B. Output: double TRUNCATED_NORMAL_AB_MOMENT, the moment of the PDF. */ { double a_h; double a_cdf; double a_pdf; double b_h; double b_cdf; double b_pdf; double ir; double irm1; double irm2; double moment; int r; if ( order < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" ); fprintf ( stderr, " ORDER < 0.\n" ); exit ( 1 ); } if ( sigma <= 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" ); fprintf ( stderr, " SIGMA <= 0.\n" ); exit ( 1 ); } if ( b <= a ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" ); fprintf ( stderr, " B <= A.\n" ); exit ( 1 ); } a_h = ( a - mu ) / sigma; a_pdf = normal_01_pdf ( a_h ); a_cdf = normal_01_cdf ( a_h ); if ( a_cdf == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" ); fprintf ( stderr, " PDF/CDF ratio fails, A_CDF too small.\n" ); fprintf ( stderr, " A_PDF = %g\n", a_pdf ); fprintf ( stderr, " A_CDF = %g\n", a_cdf ); exit ( 1 ); } b_h = ( b - mu ) / sigma; b_pdf = normal_01_pdf ( b_h ); b_cdf = normal_01_cdf ( b_h ); if ( b_cdf == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n" ); fprintf ( stderr, " PDF/CDF ratio fails, B_CDF too small.\n" ); fprintf ( stderr, " B_PDF = %g\n", b_pdf ); fprintf ( stderr, " B_CDF = %g\n", b_cdf ); exit ( 1 ); } moment = 0.0; irm2 = 0.0; irm1 = 0.0; for ( r = 0; r <= order; r++ ) { if ( r == 0 ) { ir = 1.0; } else if ( r == 1 ) { ir = - ( b_pdf - a_pdf ) / ( b_cdf - a_cdf ); } else { ir = ( double ) ( r - 1 ) * irm2 - ( pow ( b_h, r - 1 ) * b_pdf - pow ( a_h, r - 1 ) * a_pdf ) / ( b_cdf - a_cdf ); } moment = moment + r8_choose ( order, r ) * pow ( mu, order - r ) * pow ( sigma, r ) * ir; irm2 = irm1; irm1 = ir; } return moment; } /******************************************************************************/ double truncated_normal_ab_pdf ( double x, double mu, double sigma, double a, double b ) /******************************************************************************/ /* Purpose: truncated_normal_ab_pdf() evaluates the truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 14 August 2013 Author: John Burkardt Input: double X, the argument of the PDF. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, B, the lower and upper truncation limits. Output: double TRUNCATED_NORMAL_AB_PDF, the value of the PDF. */ { double alpha; double alpha_cdf; double beta; double beta_cdf; double pdf; double xi; double xi_pdf; if ( x < a ) { pdf = 0.0; } else if ( x <= b ) { alpha = ( a - mu ) / sigma; beta = ( b - mu ) / sigma; xi = ( x - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); beta_cdf = normal_01_cdf ( beta ); xi_pdf = normal_01_pdf ( xi ); pdf = xi_pdf / ( beta_cdf - alpha_cdf ) / sigma; } else { pdf = 0.0; } return pdf; } /******************************************************************************/ void truncated_normal_ab_pdf_values ( int *n_data, double *mu, double *sigma, double *a, double *b, double *x, double *fx ) /******************************************************************************/ /* Purpose: truncated_normal_ab_pdf_values(): values of the Truncated Normal PDF. Discussion: The Normal distribution, with mean Mu and standard deviation Sigma, is truncated to the interval [A,B]. Licensing: This code is distributed under the MIT license. Modified: 13 September 2013 Author: John Burkardt Reference: Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *MU, the mean of the distribution. double *SIGMA, the standard deviation of the distribution. double *A, *B, the lower and upper truncation limits. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 11 static double a_vec[N_MAX] = { 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0 }; static double b_vec[N_MAX] = { 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0 }; static double fx_vec[N_MAX] = { 0.01543301171801836, 0.01588394472270638, 0.01624375997031919, 0.01650575046469259, 0.01666496869385951, 0.01671838200940538, 0.01666496869385951, 0.01650575046469259, 0.01624375997031919, 0.01588394472270638, 0.01543301171801836 }; static double mu_vec[N_MAX] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 }; static double sigma_vec[N_MAX] = { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 }; static double x_vec[N_MAX] = { 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0, 110.0 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *b = 0.0; *mu = 0.0; *sigma = 0.0; *x = 0.0; *fx = 0.0; } else { *a = a_vec[*n_data-1]; *b = b_vec[*n_data-1]; *mu = mu_vec[*n_data-1]; *sigma = sigma_vec[*n_data-1]; *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double truncated_normal_ab_sample ( double mu, double sigma, double a, double b, int *seed ) /******************************************************************************/ /* Purpose: truncated_normal_ab_sample() samples the truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 14 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, B, the lower and upper truncation limits. int *SEED, a seed for the random number generator. Output: int *SEED: an updated seed. double TRUNCATED_NORMAL_AB_SAMPLE, a sample of the PDF. */ { double alpha; double alpha_cdf; double beta; double beta_cdf; double u; double x; double xi; double xi_cdf; alpha = ( a - mu ) / sigma; beta = ( b - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); beta_cdf = normal_01_cdf ( beta ); u = r8_uniform_01 ( seed ); xi_cdf = alpha_cdf + u * ( beta_cdf - alpha_cdf ); xi = normal_01_cdf_inv ( xi_cdf ); x = mu + sigma * xi; return x; } /******************************************************************************/ double truncated_normal_ab_variance ( double mu, double sigma, double a, double b ) /******************************************************************************/ /* Purpose: truncated_normal_ab_variance() returns the variance of the truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 14 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, B, the lower and upper truncation limits. Output: double TRUNCATED_NORMAL_AB_VARIANCE, the variance of the PDF. */ { double alpha; double alpha_cdf; double alpha_pdf; double beta; double beta_cdf; double beta_pdf; double variance; alpha = ( a - mu ) / sigma; beta = ( b - mu ) / sigma; alpha_pdf = normal_01_pdf ( alpha ); beta_pdf = normal_01_pdf ( beta ); alpha_cdf = normal_01_cdf ( alpha ); beta_cdf = normal_01_cdf ( beta ); variance = sigma * sigma * ( 1.0 + ( alpha * alpha_pdf - beta * beta_pdf ) / ( beta_cdf - alpha_cdf ) - pow ( ( alpha_pdf - beta_pdf ) / ( beta_cdf - alpha_cdf ), 2 ) ); return variance; } /******************************************************************************/ double truncated_normal_a_cdf ( double x, double mu, double sigma, double a ) /******************************************************************************/ /* Purpose: truncated_normal_a_cdf() evaluates the lower truncated Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 24 January 2017 Author: John Burkardt Input: double X, the argument of the CDF. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, the lower truncation limit. Output: double TRUNCATED_NORMAL_A_CDF, the value of the CDF. */ { double alpha; double alpha_cdf; double cdf; double xi; double xi_cdf; if ( x < a ) { cdf = 0.0; } else { alpha = ( a - mu ) / sigma; xi = ( x - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); xi_cdf = normal_01_cdf ( xi ); cdf = ( xi_cdf - alpha_cdf ) / ( 1.0 - alpha_cdf ); } return cdf; } /******************************************************************************/ void truncated_normal_a_cdf_values ( int *n_data, double *mu, double *sigma, double *a, double *x, double *fx ) /******************************************************************************/ /* Purpose: truncated_normal_a_cdf_values(): values of the Lower Truncated Normal CDF. Discussion: The Normal distribution, with mean Mu and standard deviation Sigma, is truncated to the interval [A,+oo). Licensing: This code is distributed under the MIT license. Modified: 14 September 2013 Author: John Burkardt Reference: Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *MU, the mean of the distribution. double *SIGMA, the standard deviation of the distribution. double *A, the lower truncation limit. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 11 static double a_vec[N_MAX] = { 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0 }; static double fx_vec[N_MAX] = { 0.3293202045481688, 0.3599223134505957, 0.3913175216041539, 0.4233210140873113, 0.4557365629792204, 0.4883601253415709, 0.5209836877039214, 0.5533992365958304, 0.5854027290789878, 0.6167979372325460, 0.6474000461349729 }; static double mu_vec[N_MAX] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 }; static double sigma_vec[N_MAX] = { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 }; static double x_vec[N_MAX] = { 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0, 110.0 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *mu = 0.0; *sigma = 0.0; *x = 0.0; *fx = 0.0; } else { *a = a_vec[*n_data-1]; *mu = mu_vec[*n_data-1]; *sigma = sigma_vec[*n_data-1]; *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double truncated_normal_a_cdf_inv ( double cdf, double mu, double sigma, double a ) /******************************************************************************/ /* Purpose: truncated_normal_a_cdf_inv() inverts the lower truncated Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double CDF, the value of the CDF. 0.0D+00 <= CDF <= 1.0. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, the lower truncation limit. Output: double TRUNCATED_NORMAL_A_CDF_INV, the corresponding argument. */ { double alpha; double alpha_cdf; double x; double xi; double xi_cdf; if ( cdf < 0.0 || 1.0 < cdf ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_A_CDF_INV - Fatal error!\n" ); fprintf ( stderr, " CDF < 0 or 1 < CDF.\n" ); exit ( 1 ); } alpha = ( a - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); xi_cdf = ( 1.0 - alpha_cdf ) * cdf + alpha_cdf; xi = normal_01_cdf_inv ( xi_cdf ); x = mu + sigma * xi; return x; } /******************************************************************************/ double truncated_normal_a_mean ( double mu, double sigma, double a ) /******************************************************************************/ /* Purpose: truncated_normal_a_mean() returns the mean of the lower truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, the lower truncation limit. Output: double TRUNCATED_NORMAL_A_MEAN, the mean of the PDF. */ { double alpha; double alpha_cdf; double alpha_pdf; double mean; alpha = ( a - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); alpha_pdf = normal_01_pdf ( alpha ); mean = mu + sigma * alpha_pdf / ( 1.0 - alpha_cdf ); return mean; } /******************************************************************************/ double truncated_normal_a_moment ( int order, double mu, double sigma, double a ) /******************************************************************************/ /* Purpose: truncated_normal_a_moment(): moments of the lower truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 11 September 2013 Author: John Burkardt Reference: Phoebus Dhrymes, Moments of Truncated Normal Distributions, May 2005. Input: int ORDER, the order of the moment. 0 <= ORDER. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, the lower truncation limit. Output: double TRUNCATED_NORMAL_A_MOMENT, the moment of the PDF. */ { double moment; moment = r8_mop ( order ) * truncated_normal_b_moment ( order, - mu, sigma, - a ); return moment; } /******************************************************************************/ double truncated_normal_a_pdf ( double x, double mu, double sigma, double a ) /******************************************************************************/ /* Purpose: truncated_normal_a_pdf() evaluates the lower truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 24 January 2017 Author: John Burkardt Input: double X, the argument of the PDF. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, the lower truncation limit. Output: double TRUNCATED_NORMAL_A_PDF, the value of the PDF. */ { double alpha; double alpha_cdf; double pdf; double xi; double xi_pdf; if ( x < a ) { pdf = 0.0; } else { alpha = ( a - mu ) / sigma; xi = ( x - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); xi_pdf = normal_01_pdf ( xi ); pdf = xi_pdf / ( 1.0 - alpha_cdf ) / sigma; } return pdf; } /******************************************************************************/ void truncated_normal_a_pdf_values ( int *n_data, double *mu, double *sigma, double *a, double *x, double *fx ) /******************************************************************************/ /* Purpose: truncated_normal_a_pdf_values(): values of the Lower Truncated Normal PDF. Discussion: The Normal distribution, with mean Mu and standard deviation Sigma, is truncated to the interval [A,+oo). Licensing: This code is distributed under the MIT license. Modified: 14 September 2013 Author: John Burkardt Reference: Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *MU, the mean of the distribution. double *SIGMA, the standard deviation of the distribution. double *A, the lower truncation limit. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 11 static double a_vec[N_MAX] = { 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0, 50.0 }; static double fx_vec[N_MAX] = { 0.01507373507401876, 0.01551417047139894, 0.01586560931024694, 0.01612150073158793, 0.01627701240029317, 0.01632918226724295, 0.01627701240029317, 0.01612150073158793, 0.01586560931024694, 0.01551417047139894, 0.01507373507401876 }; static double mu_vec[N_MAX] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 }; static double sigma_vec[N_MAX] = { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 }; static double x_vec[N_MAX] = { 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0, 110.0 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *a = 0.0; *mu = 0.0; *sigma = 0.0; *x = 0.0; *fx = 0.0; } else { *a = a_vec[*n_data-1]; *mu = mu_vec[*n_data-1]; *sigma = sigma_vec[*n_data-1]; *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double truncated_normal_a_sample ( double mu, double sigma, double a, int *seed ) /******************************************************************************/ /* Purpose: truncated_normal_a_sample() samples the lower truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, the lower truncation limit. int *SEED, a seed for the random number generator. Output: int *SEED: an updated seed. double TRUNCATED_NORMAL_A_SAMPLE, a sample of the PDF. */ { double alpha; double alpha_cdf; double u; double x; double xi; double xi_cdf; alpha = ( a - mu ) / sigma; alpha_cdf = normal_01_cdf ( alpha ); u = r8_uniform_01 ( seed ); xi_cdf = alpha_cdf + u * ( 1.0 - alpha_cdf ); xi = normal_01_cdf_inv ( xi_cdf ); x = mu + sigma * xi; return x; } /******************************************************************************/ double truncated_normal_a_variance ( double mu, double sigma, double a ) /******************************************************************************/ /* Purpose: truncated_normal_a_variance(): variance of the lower truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double A, the lower truncation limit. Output: double TRUNCATED_NORMAL_A_VARIANCE, the variance of the PDF. */ { double alpha; double alpha_cdf; double alpha_pdf; double variance; alpha = ( a - mu ) / sigma; alpha_pdf = normal_01_pdf ( alpha ); alpha_cdf = normal_01_cdf ( alpha ); variance = sigma * sigma * ( 1.0 + alpha * alpha_pdf / ( 1.0 - alpha_cdf ) - pow ( alpha_pdf / ( 1.0 - alpha_cdf ), 2 ) ); return variance; } /******************************************************************************/ double truncated_normal_b_cdf ( double x, double mu, double sigma, double b ) /******************************************************************************/ /* Purpose: truncated_normal_b_cdf() evaluates the upper truncated Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 24 January 2017 Author: John Burkardt Input: double X, the argument of the CDF. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double B, the upper truncation limit. Output: double TRUNCATED_NORMAL_B_CDF, the value of the CDF. */ { double beta; double beta_cdf; double cdf; double xi; double xi_cdf; if ( x <= b ) { beta = ( b - mu ) / sigma; xi = ( x - mu ) / sigma; beta_cdf = normal_01_cdf ( beta ); xi_cdf = normal_01_cdf ( xi ); cdf = xi_cdf / beta_cdf; } else { cdf = 1.0; } return cdf; } /******************************************************************************/ void truncated_normal_b_cdf_values ( int *n_data, double *mu, double *sigma, double *b, double *x, double *fx ) /******************************************************************************/ /* Purpose: truncated_normal_b_cdf_values(): values of the Upper Truncated Normal CDF. Discussion: The Normal distribution, with mean Mu and standard deviation Sigma, is truncated to the interval [A,B]. Licensing: This code is distributed under the MIT license. Modified: 14 September 2013 Author: John Burkardt Reference: Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *MU, the mean of the distribution. double *SIGMA, the standard deviation of the distribution. double *B, the upper truncation limit. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 11 static double b_vec[N_MAX] = { 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0 }; static double fx_vec[N_MAX] = { 0.3525999538650271, 0.3832020627674540, 0.4145972709210122, 0.4466007634041696, 0.4790163122960786, 0.5116398746584291, 0.5442634370207796, 0.5766789859126887, 0.6086824783958461, 0.6400776865494043, 0.6706797954518312 }; static double mu_vec[N_MAX] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 }; static double sigma_vec[N_MAX] = { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 }; static double x_vec[N_MAX] = { 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0, 110.0 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *b = 0.0; *mu = 0.0; *sigma = 0.0; *x = 0.0; *fx = 0.0; } else { *b = b_vec[*n_data-1]; *mu = mu_vec[*n_data-1]; *sigma = sigma_vec[*n_data-1]; *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double truncated_normal_b_cdf_inv ( double cdf, double mu, double sigma, double b ) /******************************************************************************/ /* Purpose: truncated_normal_b_cdf_inv() inverts the upper truncated Normal CDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double CDF, the value of the CDF. 0.0D+00 <= CDF <= 1.0. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double B, the upper truncation limit. Output: double TRUNCATED_NORMAL_B_CDF_INV, the corresponding argument. */ { double beta; double beta_cdf; double x; double xi; double xi_cdf; if ( cdf < 0.0 || 1.0 < cdf ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_B_CDF_INV - Fatal error!\n" ); fprintf ( stderr, " CDF < 0 or 1 < CDF.\n" ); exit ( 1 ); } beta = ( b - mu ) / sigma; beta_cdf = normal_01_cdf ( beta ); xi_cdf = beta_cdf * cdf; xi = normal_01_cdf_inv ( xi_cdf ); x = mu + sigma * xi; return x; } /******************************************************************************/ double truncated_normal_b_mean ( double mu, double sigma, double b ) /******************************************************************************/ /* Purpose: truncated_normal_b_mean() returns the mean of the upper truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double B, the upper truncation limit. Output: double TRUNCATED_NORMAL_B_MEAN, the mean of the PDF. */ { double beta; double beta_cdf; double beta_pdf; double mean; beta = ( b - mu ) / sigma; beta_cdf = normal_01_cdf ( beta ); beta_pdf = normal_01_pdf ( beta ); mean = mu - sigma * beta_pdf / beta_cdf; return mean; } /******************************************************************************/ double truncated_normal_b_moment ( int order, double mu, double sigma, double b ) /******************************************************************************/ /* Purpose: truncated_normal_b_moment(): moments of the upper truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 11 September 2013 Author: John Burkardt Reference: Phoebus Dhrymes, Moments of Truncated Normal Distributions, May 2005. Input: int ORDER, the order of the moment. 0 <= ORDER. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double B, the upper truncation limit. Output: double TRUNCATED_NORMAL_B_MOMENT, the moment of the PDF. */ { double f; double h; double h_cdf; double h_pdf; double ir; double irm1; double irm2; double moment; int r; if ( order < 0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n" ); fprintf ( stderr, " ORDER < 0.\n" ); exit ( 1 ); } h = ( b - mu ) / sigma; h_pdf = normal_01_pdf ( h ); h_cdf = normal_01_cdf ( h ); if ( h_cdf == 0.0 ) { fprintf ( stderr, "\n" ); fprintf ( stderr, "TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n" ); fprintf ( stderr, " CDF((B-MU)/S) = 0.\n" ); exit ( 1 ); } f = h_pdf / h_cdf; moment = 0.0; irm2 = 0.0; irm1 = 0.0; for ( r = 0; r <= order; r++ ) { if ( r == 0 ) { ir = 1.0; } else if ( r == 1 ) { ir = - f; } else { ir = - pow ( h, r - 1 ) * f + ( double ) ( r - 1 ) * irm2; } moment = moment + r8_choose ( order, r ) * pow ( mu, order - r ) * pow ( sigma, r ) * ir; irm2 = irm1; irm1 = ir; } return moment; } /******************************************************************************/ double truncated_normal_b_pdf ( double x, double mu, double sigma, double b ) /******************************************************************************/ /* Purpose: truncated_normal_b_pdf() evaluates the upper truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 24 January 2017 Author: John Burkardt Input: double X, the argument of the PDF. double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double B, the upper truncation limit. Output: double TRUNCATED_NORMAL_B_PDF, the value of the PDF. */ { double beta; double beta_cdf; double pdf; double xi; double xi_pdf; if ( x <= b ) { beta = ( b - mu ) / sigma; xi = ( x - mu ) / sigma; beta_cdf = normal_01_cdf ( beta ); xi_pdf = normal_01_pdf ( xi ); pdf = xi_pdf / beta_cdf / sigma; } else { pdf = 0.0; } return pdf; } /******************************************************************************/ void truncated_normal_b_pdf_values ( int *n_data, double *mu, double *sigma, double *b, double *x, double *fx ) /******************************************************************************/ /* Purpose: truncated_normal_b_pdf_values(): values of the Upper Truncated Normal PDF. Discussion: The Normal distribution, with mean Mu and standard deviation Sigma, is truncated to the interval (-oo,B]. Licensing: This code is distributed under the MIT license. Modified: 14 September 2013 Author: John Burkardt Reference: Stephen Wolfram, The Mathematica Book, Fourth Edition, Cambridge University Press, 1999, ISBN: 0-521-64314-7, LC: QA76.95.W65. Input: int *N_DATA: The user sets N_DATA to 0 before the first call. Output: int *N_DATA: On each call, the routine increments N_DATA by 1, and returns the corresponding data; when there is no more data, the output value of N_DATA will be 0 again. double *MU, the mean of the distribution. double *SIGMA, the standard deviation of the distribution. double *B, the upper truncation limit. double *X, the argument of the function. double *FX, the value of the function. */ { # define N_MAX 11 static double b_vec[N_MAX] = { 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0, 150.0 }; static double fx_vec[N_MAX] = { 0.01507373507401876, 0.01551417047139894, 0.01586560931024694, 0.01612150073158793, 0.01627701240029317, 0.01632918226724295, 0.01627701240029317, 0.01612150073158793, 0.01586560931024694, 0.01551417047139894, 0.01507373507401876 }; static double mu_vec[N_MAX] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 }; static double sigma_vec[N_MAX] = { 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0, 25.0 }; static double x_vec[N_MAX] = { 90.0, 92.0, 94.0, 96.0, 98.0, 100.0, 102.0, 104.0, 106.0, 108.0, 110.0 }; if ( *n_data < 0 ) { *n_data = 0; } *n_data = *n_data + 1; if ( N_MAX < *n_data ) { *n_data = 0; *b = 0.0; *mu = 0.0; *sigma = 0.0; *x = 0.0; *fx = 0.0; } else { *b = b_vec[*n_data-1]; *mu = mu_vec[*n_data-1]; *sigma = sigma_vec[*n_data-1]; *x = x_vec[*n_data-1]; *fx = fx_vec[*n_data-1]; } return; # undef N_MAX } /******************************************************************************/ double truncated_normal_b_sample ( double mu, double sigma, double b, int *seed ) /******************************************************************************/ /* Purpose: truncated_normal_b_sample() samples the upper truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double B, the upper truncation limit. int *SEED, a seed for the random number generator. Output: int *SEED: an updated seed. double TRUNCATED_NORMAL_B_SAMPLE, a sample of the PDF. */ { double beta; double beta_cdf; double u; double x; double xi; double xi_cdf; beta = ( b - mu ) / sigma; beta_cdf = normal_01_cdf ( beta ); u = r8_uniform_01 ( seed ); xi_cdf = u * beta_cdf; xi = normal_01_cdf_inv ( xi_cdf ); x = mu + sigma * xi; return x; } /******************************************************************************/ double truncated_normal_b_variance ( double mu, double sigma, double b ) /******************************************************************************/ /* Purpose: truncated_normal_b_variance(): variance of the upper truncated Normal PDF. Licensing: This code is distributed under the MIT license. Modified: 21 August 2013 Author: John Burkardt Input: double MU, SIGMA, the mean and standard deviation of the parent Normal distribution. double B, the upper truncation limit. Output: double TRUNCATED_NORMAL_B_VARIANCE, the variance of the PDF. */ { double beta; double beta_cdf; double beta_pdf; double variance; beta = ( b - mu ) / sigma; beta_pdf = normal_01_pdf ( beta ); beta_cdf = normal_01_cdf ( beta ); variance = sigma * sigma * ( 1.0 - beta * beta_pdf / beta_cdf - pow ( beta_pdf / beta_cdf, 2 ) ); return variance; }