|
|
@@ -0,0 +1,3985 @@
|
|
|
+#include <cmath>
|
|
|
+#include <cstdlib>
|
|
|
+#include <ctime>
|
|
|
+#include <iomanip>
|
|
|
+#include <iostream>
|
|
|
+#include <string>
|
|
|
+
|
|
|
+using namespace std;
|
|
|
+
|
|
|
+#include "truncated_normal.hpp"
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+int i4_uniform_ab(int a, int b, int &seed)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 02 October 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int A, B, the limits of the interval.
|
|
|
+//
|
|
|
+// Input/output, int &SEED, the "seed" value, which should NOT be 0.
|
|
|
+// On output, SEED has been updated.
|
|
|
+//
|
|
|
+// Output, 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "I4_UNIFORM_AB - Fatal error!\n";
|
|
|
+ cerr << " Input value of SEED = 0.\n";
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+ //
|
|
|
+ // Guarantee 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);
|
|
|
+ //
|
|
|
+ // Use rounding to convert R to an integer between A and B.
|
|
|
+ //
|
|
|
+ value = round(r);
|
|
|
+ //
|
|
|
+ // Guarantee A <= VALUE <= B.
|
|
|
+ //
|
|
|
+ if (value < a) {
|
|
|
+ value = a;
|
|
|
+ }
|
|
|
+ if (b < value) {
|
|
|
+ value = b;
|
|
|
+ }
|
|
|
+
|
|
|
+ return value;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_01_cdf(double x)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_01_CDF evaluates the Normal 01 CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_01_cdf_inv(double p)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 27 December 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// Original FORTRAN77 version by Michael Wichura.
|
|
|
+// C++ 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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+2,
|
|
|
+ 1.9715909503065514427E+3, 1.3731693765509461125E+4,
|
|
|
+ 4.5921953931549871457E+4, 6.7265770927008700853E+4,
|
|
|
+ 3.3430575583588128105E+4, 2.5090809287301226727E+3};
|
|
|
+ double b[8] = {1.0,
|
|
|
+ 4.2313330701600911252E+1,
|
|
|
+ 6.8718700749205790830E+2,
|
|
|
+ 5.3941960214247511077E+3,
|
|
|
+ 2.1213794301586595867E+4,
|
|
|
+ 3.9307895800092710610E+4,
|
|
|
+ 2.8729085735721942674E+4,
|
|
|
+ 5.2264952788528545610E+3};
|
|
|
+ double c[8] = {1.42343711074968357734, 4.63033784615654529590,
|
|
|
+ 5.76949722146069140550, 3.64784832476320460504,
|
|
|
+ 1.27045825245236838258, 2.41780725177450611770E-1,
|
|
|
+ 2.27238449892691845833E-2, 7.74545014278341407640E-4};
|
|
|
+ double const1 = 0.180625;
|
|
|
+ double const2 = 1.6;
|
|
|
+ double d[8] = {1.0,
|
|
|
+ 2.05319162663775882187,
|
|
|
+ 1.67638483018380384940,
|
|
|
+ 6.89767334985100004550E-1,
|
|
|
+ 1.48103976427480074590E-1,
|
|
|
+ 1.51986665636164571966E-2,
|
|
|
+ 5.47593808499534494600E-4,
|
|
|
+ 1.05075007164441684324E-9};
|
|
|
+ double e[8] = {6.65790464350110377720, 5.46378491116411436990,
|
|
|
+ 1.78482653991729133580, 2.96560571828504891230E-1,
|
|
|
+ 2.65321895265761230930E-2, 1.24266094738807843860E-3,
|
|
|
+ 2.71155556874348757815E-5, 2.01033439929228813265E-7};
|
|
|
+ double f[8] = {1.0,
|
|
|
+ 5.99832206555887937690E-1,
|
|
|
+ 1.36929880922735805310E-1,
|
|
|
+ 1.48753612908506148525E-2,
|
|
|
+ 7.86869131145613259100E-4,
|
|
|
+ 1.84631831751005468180E-5,
|
|
|
+ 1.42151175831644588870E-7,
|
|
|
+ 2.04426310338993978564E-15};
|
|
|
+ double q;
|
|
|
+ double r;
|
|
|
+ double split1 = 0.425;
|
|
|
+ double split2 = 5.0;
|
|
|
+ double value;
|
|
|
+
|
|
|
+ if (p <= 0.0) {
|
|
|
+ value = -r8_huge();
|
|
|
+ return value;
|
|
|
+ }
|
|
|
+
|
|
|
+ if (1.0 <= p) {
|
|
|
+ value = r8_huge();
|
|
|
+ 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 = r8_huge();
|
|
|
+ } 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void normal_01_cdf_values(int &n_data, double &x, double &fx)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, double &X, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_01_mean()
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_01_MEAN returns the mean of the Normal 01 PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 18 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Output, double MEAN, the mean of the PDF.
|
|
|
+//
|
|
|
+{
|
|
|
+ double mean;
|
|
|
+
|
|
|
+ mean = 0.0;
|
|
|
+
|
|
|
+ return mean;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_01_moment(int order)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_01_MOMENT evaluates moments of the Normal 01 PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 31 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_01_pdf(double x)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 18 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the PDF.
|
|
|
+//
|
|
|
+// Output, double 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_01_sample(int &seed)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 26 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &SEED, a seed for the random number generator.
|
|
|
+//
|
|
|
+// Output, double NORMAL_01_SAMPLE, a normally distributed random value.
|
|
|
+//
|
|
|
+{
|
|
|
+ double r1;
|
|
|
+ double r2;
|
|
|
+ const double r8_pi = 3.14159265358979323;
|
|
|
+ double x;
|
|
|
+
|
|
|
+ r1 = r8_uniform_01(seed);
|
|
|
+ r2 = r8_uniform_01(seed);
|
|
|
+ x = sqrt(-2.0 * log(r1)) * cos(2.0 * r8_pi * r2);
|
|
|
+
|
|
|
+ return x;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_01_variance()
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_01_VARIANCE returns the variance of the Normal 01 PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 18 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Output, double VARIANCE, the variance of the PDF.
|
|
|
+//
|
|
|
+{
|
|
|
+ double variance;
|
|
|
+
|
|
|
+ variance = 1.0;
|
|
|
+
|
|
|
+ return variance;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_cdf(double x, double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_MS_CDF evaluates the Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 19 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the CDF.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_cdf_inv(double cdf, double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_MS_CDF_INV inverts the Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 19 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Reference:
|
|
|
+//
|
|
|
+// Algorithm AS 111,
|
|
|
+// Applied Statistics,
|
|
|
+// Volume 26, pages 118-121, 1977.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double CDF, the value of the CDF.
|
|
|
+// 0.0 <= CDF <= 1.0.
|
|
|
+//
|
|
|
+// Input, 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) {
|
|
|
+ cout << "\n";
|
|
|
+ cout << "NORMAL_MS_CDF_INV - Fatal error!\n";
|
|
|
+ cout << " CDF < 0 or 1 < CDF.\n";
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
+ x2 = normal_01_cdf_inv(cdf);
|
|
|
+
|
|
|
+ x = mu + sigma * x2;
|
|
|
+
|
|
|
+ return x;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_mean(double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_MS_MEAN returns the mean of the Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 19 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_moment(int order, double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 31 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int ORDER, the order of the moment.
|
|
|
+// 0 <= ORDER.
|
|
|
+//
|
|
|
+// Input, double MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_moment_central(int order, double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_MS_MOMENT_CENTRAL evaluates central moments of the Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 31 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int ORDER, the order of the moment.
|
|
|
+// 0 <= ORDER.
|
|
|
+//
|
|
|
+// Input, double MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_moment_central_values(int order, double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 31 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int ORDER, the order of the moment.
|
|
|
+// 0 <= ORDER <= 10.
|
|
|
+//
|
|
|
+// Input, double MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Input, 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 {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "NORMAL_MS_MOMENT_CENTRAL_VALUES - Fatal error!\n";
|
|
|
+ cerr << " Only ORDERS 0 through 10 are available.\n";
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
+ return value;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_moment_values(int order, double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 31 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int ORDER, the order of the moment.
|
|
|
+// 0 <= ORDER <= 8.
|
|
|
+//
|
|
|
+// Input, double MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Input, 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 {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "NORMAL_MS_MOMENT_VALUES - Fatal error!\n";
|
|
|
+ cerr << " Only ORDERS 0 through 8 are available.\n";
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
+ return value;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_pdf(double x, double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 19 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the PDF.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_sample(double mu, double sigma, int &seed)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_MS_SAMPLE samples the Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 20 November 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the parameters of the PDF.
|
|
|
+// 0.0 < SIGMA.
|
|
|
+//
|
|
|
+// Input/output, int &SEED, a seed for the random number generator.
|
|
|
+//
|
|
|
+// Output, double NORMAL_MS_SAMPLE, a sample of the PDF.
|
|
|
+//
|
|
|
+{
|
|
|
+ double x;
|
|
|
+
|
|
|
+ x = normal_01_sample(seed);
|
|
|
+
|
|
|
+ x = mu + sigma * x;
|
|
|
+
|
|
|
+ return x;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double normal_ms_variance(double mu, double sigma)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// NORMAL_MS_VARIANCE returns the variance of the Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 19 September 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8_abs(double x)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8_ABS returns the absolute value of an R8.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 November 2006
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the quantity whose absolute value is desired.
|
|
|
+//
|
|
|
+// Output, double R8_ABS, the absolute value of X.
|
|
|
+//
|
|
|
+{
|
|
|
+ double value;
|
|
|
+
|
|
|
+ if (0.0 <= x) {
|
|
|
+ value = x;
|
|
|
+ } else {
|
|
|
+ value = -x;
|
|
|
+ }
|
|
|
+ return value;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8_choose(int n, int k)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 09 June 2013
|
|
|
+//
|
|
|
+// 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8_factorial2(int n)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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)
|
|
|
+//
|
|
|
+// Example:
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 22 January 2008
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void r8_factorial2_values(int &n_data, int &n, double &f)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, int &N, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8_huge()
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8_HUGE returns a "huge" R8.
|
|
|
+//
|
|
|
+// Discussion:
|
|
|
+//
|
|
|
+// HUGE_VAL is the largest representable legal real number, and is usually
|
|
|
+// defined in math.h, or sometimes in stdlib.h.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 08 May 2003
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Output, double R8_HUGE, a "huge" real value.
|
|
|
+//
|
|
|
+{
|
|
|
+ return HUGE_VAL;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8_log_2(double x)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8_LOG_2 returns the logarithm base 2 of the absolute value of an R8.
|
|
|
+//
|
|
|
+// Discussion:
|
|
|
+//
|
|
|
+// value = Log ( |X| ) / Log ( 2.0 )
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 22 March 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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 = -r8_huge();
|
|
|
+ } else {
|
|
|
+ value = log(fabs(x)) / log(2.0);
|
|
|
+ }
|
|
|
+
|
|
|
+ return value;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8_mop(int i)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 16 November 2007
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8_uniform_01(int &seed)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &SEED, the "seed" value. Normally, this
|
|
|
+// value should not be 0. On output, SEED has been updated.
|
|
|
+//
|
|
|
+// Output, 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;
|
|
|
+ }
|
|
|
+
|
|
|
+ r = (double)(seed)*4.656612875E-10;
|
|
|
+
|
|
|
+ return r;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void r8poly_print(int n, double a[], string title)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8POLY_PRINT prints out a polynomial.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 26 February 2015
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int N, the dimension of A.
|
|
|
+//
|
|
|
+// Input, double A[N+1], the polynomial coefficients.
|
|
|
+// A(0) is the constant term and
|
|
|
+// A(N) is the coefficient of X^N.
|
|
|
+//
|
|
|
+// Input, string TITLE, a title.
|
|
|
+//
|
|
|
+{
|
|
|
+ int i;
|
|
|
+ double mag;
|
|
|
+ int n2;
|
|
|
+ char plus_minus;
|
|
|
+
|
|
|
+ cout << "\n";
|
|
|
+ cout << title << "\n";
|
|
|
+ cout << "\n";
|
|
|
+
|
|
|
+ if (n <= 0) {
|
|
|
+ cout << " p(x) = 0\n";
|
|
|
+ return;
|
|
|
+ }
|
|
|
+
|
|
|
+ if (a[n] < 0.0) {
|
|
|
+ plus_minus = '-';
|
|
|
+ } else {
|
|
|
+ plus_minus = ' ';
|
|
|
+ }
|
|
|
+
|
|
|
+ mag = fabs(a[n]);
|
|
|
+
|
|
|
+ if (2 <= n) {
|
|
|
+ cout << " p(x) = " << plus_minus << setw(14) << mag << " * x ^ " << n
|
|
|
+ << "\n";
|
|
|
+ } else if (n == 1) {
|
|
|
+ cout << " p(x) = " << plus_minus << setw(14) << mag << " * x\n";
|
|
|
+ } else if (n == 0) {
|
|
|
+ cout << " p(x) = " << plus_minus << setw(14) << mag << "\n";
|
|
|
+ }
|
|
|
+
|
|
|
+ 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) {
|
|
|
+ cout << " " << plus_minus << setw(14) << mag << " * x ^ " << i
|
|
|
+ << "\n";
|
|
|
+ } else if (i == 1) {
|
|
|
+ cout << " " << plus_minus << setw(14) << mag << " * x\n";
|
|
|
+ } else if (i == 0) {
|
|
|
+ cout << " " << plus_minus << setw(14) << mag << "\n";
|
|
|
+ }
|
|
|
+ }
|
|
|
+ }
|
|
|
+
|
|
|
+ return;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8poly_value_horner(int m, double c[], double x)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 02 January 2015
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int M, the degree of the polynomial.
|
|
|
+//
|
|
|
+// Input, double C[M+1], the coefficients of the polynomial.
|
|
|
+// A[0] is the constant term.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double *r8vec_linspace_new(int n, double a_first, double a_last)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 29 March 2011
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int N, the number of entries in the vector.
|
|
|
+//
|
|
|
+// Input, double A_FIRST, A_LAST, the first and last entries.
|
|
|
+//
|
|
|
+// Output, double R8VEC_LINSPACE_NEW[N], a vector of linearly spaced data.
|
|
|
+//
|
|
|
+{
|
|
|
+ double *a;
|
|
|
+ int i;
|
|
|
+
|
|
|
+ a = new double[n];
|
|
|
+
|
|
|
+ if (n == 1) {
|
|
|
+ a[0] = (a_first + a_last) / 2.0;
|
|
|
+ } else {
|
|
|
+ for (i = 0; i < n; i++) {
|
|
|
+ a[i] = ((double)(n - 1 - i) * a_first + (double)(i)*a_last) /
|
|
|
+ (double)(n - 1);
|
|
|
+ }
|
|
|
+ }
|
|
|
+ return a;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8vec_max(int n, double *dvec)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8VEC_MAX returns the value of the maximum element in an R8VEC.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 15 October 1998
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int N, the number of entries in the array.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8vec_mean(int n, double x[])
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8VEC_MEAN returns the mean of an R8VEC.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 01 May 1999
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int N, the number of entries in the vector.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8vec_min(int n, double *dvec)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8VEC_MIN returns the value of the minimum element in an R8VEC.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 15 October 1998
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int N, the number of entries in the array.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void r8vec_print(int n, double a[], string title)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8VEC_PRINT prints an R8VEC.
|
|
|
+//
|
|
|
+// Discussion:
|
|
|
+//
|
|
|
+// An R8VEC is a vector of R8's.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 16 August 2004
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int N, the number of components of the vector.
|
|
|
+//
|
|
|
+// Input, double A[N], the vector to be printed.
|
|
|
+//
|
|
|
+// Input, string TITLE, a title.
|
|
|
+//
|
|
|
+{
|
|
|
+ int i;
|
|
|
+
|
|
|
+ cout << "\n";
|
|
|
+ cout << title << "\n";
|
|
|
+ cout << "\n";
|
|
|
+ for (i = 0; i < n; i++) {
|
|
|
+ cout << " " << setw(8) << i << ": " << setw(14) << a[i] << "\n";
|
|
|
+ }
|
|
|
+
|
|
|
+ return;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double r8vec_variance(int n, double x[])
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// R8VEC_VARIANCE returns the variance of an R8VEC.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 01 May 1999
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int N, the number of entries in the vector.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void timestamp()
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TIMESTAMP prints the current YMDHMS date as a time stamp.
|
|
|
+//
|
|
|
+// Example:
|
|
|
+//
|
|
|
+// 31 May 2001 09:45:54 AM
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 08 July 2009
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// None
|
|
|
+//
|
|
|
+{
|
|
|
+#define TIME_SIZE 40
|
|
|
+
|
|
|
+ static char time_buffer[TIME_SIZE];
|
|
|
+ const struct std::tm *tm_ptr;
|
|
|
+ size_t len;
|
|
|
+ std::time_t now;
|
|
|
+
|
|
|
+ now = std::time(NULL);
|
|
|
+ tm_ptr = std::localtime(&now);
|
|
|
+
|
|
|
+ len = std::strftime(time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm_ptr);
|
|
|
+
|
|
|
+ std::cout << time_buffer << "\n";
|
|
|
+
|
|
|
+ return;
|
|
|
+#undef TIME_SIZE
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_ab_cdf(double x, double mu, double sigma, double a,
|
|
|
+ double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_AB_CDF evaluates the truncated Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the CDF.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+
|
|
|
+ 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);
|
|
|
+
|
|
|
+ return cdf;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void truncated_normal_ab_cdf_values(int &n_data, double &mu, double &sigma,
|
|
|
+ double &a, double &b, double &x, double &fx)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, double &MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Output, double &SIGMA, the standard deviation of the distribution.
|
|
|
+//
|
|
|
+// Output, double &A, &B, the lower and upper truncation limits.
|
|
|
+//
|
|
|
+// Output, double &X, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_ab_cdf_inv(double cdf, double mu, double sigma,
|
|
|
+ double a, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_AB_CDF_INV inverts the truncated Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double CDF, the value of the CDF.
|
|
|
+// 0.0 <= CDF <= 1.0.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_AB_CDF_INV - Fatal error!\n";
|
|
|
+ cerr << " 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_ab_mean(double mu, double sigma, double a, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_AB_MEAN returns the mean of the truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_ab_moment(int order, double mu, double sigma, double a,
|
|
|
+ double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_AB_MOMENT: moments of the truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 11 September 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Reference:
|
|
|
+//
|
|
|
+// Phoebus Dhrymes,
|
|
|
+// Moments of Truncated Normal Distributions,
|
|
|
+// May 2005.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int ORDER, the order of the moment.
|
|
|
+// 0 <= ORDER.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+// 0.0 < S.
|
|
|
+//
|
|
|
+// Input, double A, B, the lower and upper truncation limits.
|
|
|
+// A < B.
|
|
|
+//
|
|
|
+// Output, double TRUNCATED_NORMAL_AB_MOMENT, the moment of the PDF.
|
|
|
+//
|
|
|
+{
|
|
|
+ double a_cdf;
|
|
|
+ double a_h;
|
|
|
+ double a_pdf;
|
|
|
+ double b_cdf;
|
|
|
+ double b_h;
|
|
|
+ double b_pdf;
|
|
|
+ double ir;
|
|
|
+ double irm1;
|
|
|
+ double irm2;
|
|
|
+ double moment;
|
|
|
+ int r;
|
|
|
+
|
|
|
+ if (order < 0) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n";
|
|
|
+ cerr << " ORDER < 0.\n";
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
+ if (sigma <= 0.0) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n";
|
|
|
+ cerr << " SIGMA <= 0.0.\n";
|
|
|
+ exit(1);
|
|
|
+ }
|
|
|
+
|
|
|
+ if (b <= a) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n";
|
|
|
+ cerr << " 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n";
|
|
|
+ cerr << " PDF/CDF ratio fails, because A_CDF too small.\n";
|
|
|
+ cerr << " A_PDF = " << a_pdf << "\n";
|
|
|
+ cerr << " A_CDF = " << a_cdf << "\n";
|
|
|
+ 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_AB_MOMENT - Fatal error!\n";
|
|
|
+ cerr << " PDF/CDF ratio fails, because B_CDF too small.\n";
|
|
|
+ cerr << " B_PDF = " << b_pdf << "\n";
|
|
|
+ cerr << " B_CDF = " << b_cdf << "\n";
|
|
|
+ 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_ab_pdf(double x, double mu, double sigma, double a,
|
|
|
+ double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_AB_PDF evaluates the truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the PDF.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+
|
|
|
+ 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;
|
|
|
+
|
|
|
+ return pdf;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void truncated_normal_ab_pdf_values(int &n_data, double &mu, double &sigma,
|
|
|
+ double &a, double &b, double &x, double &fx)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, double &MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Output, double &SIGMA, the standard deviation of the distribution.
|
|
|
+//
|
|
|
+// Output, double &A, &B, the lower and upper truncation limits.
|
|
|
+//
|
|
|
+// Output, double &X, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_ab_sample(double mu, double sigma, double a, double b,
|
|
|
+ int &seed)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_AB_SAMPLE samples the truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, double A, B, the lower and upper truncation limits.
|
|
|
+//
|
|
|
+// Input/output, int &SEED, a seed for the random number
|
|
|
+// generator.
|
|
|
+//
|
|
|
+// Output, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_ab_variance(double mu, double sigma, double a, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_AB_VARIANCE returns the variance of the truncated Normal
|
|
|
+// PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_a_cdf(double x, double mu, double sigma, double a)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_A_CDF evaluates the lower truncated Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the CDF.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+
|
|
|
+ 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void truncated_normal_a_cdf_values(int &n_data, double &mu, double &sigma,
|
|
|
+ double &a, double &x, double &fx)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, double &MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Output, double &SIGMA, the standard deviation of the distribution.
|
|
|
+//
|
|
|
+// Output, double &A, the lower truncation limit.
|
|
|
+//
|
|
|
+// Output, double &X, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_a_cdf_inv(double cdf, double mu, double sigma, double a)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_A_CDF_INV inverts the lower truncated Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double CDF, the value of the CDF.
|
|
|
+// 0.0 <= CDF <= 1.0.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_A_CDF_INV - Fatal error!\n";
|
|
|
+ cerr << " 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_a_mean(double mu, double sigma, double a)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_A_MEAN returns the mean of the lower truncated Normal
|
|
|
+// PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_a_moment(int order, double mu, double sigma, double a)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_A_MOMENT: moments of the lower truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 11 September 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Reference:
|
|
|
+//
|
|
|
+// Phoebus Dhrymes,
|
|
|
+// Moments of Truncated Normal Distributions,
|
|
|
+// May 2005.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int ORDER, the order of the moment.
|
|
|
+// 0 <= ORDER.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_a_pdf(double x, double mu, double sigma, double a)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_A_PDF evaluates the lower truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the PDF.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+
|
|
|
+ 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void truncated_normal_a_pdf_values(int &n_data, double &mu, double &sigma,
|
|
|
+ double &a, double &x, double &fx)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, double &MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Output, double &SIGMA, the standard deviation of the distribution.
|
|
|
+//
|
|
|
+// Output, double &A, the lower truncation limit.
|
|
|
+//
|
|
|
+// Output, double &X, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_a_sample(double mu, double sigma, double a, int &seed)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_A_SAMPLE samples the lower truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, double A, the lower truncation limit.
|
|
|
+//
|
|
|
+// Input/output, int &SEED, a seed for the random number
|
|
|
+// generator.
|
|
|
+//
|
|
|
+// Output, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_a_variance(double mu, double sigma, double a)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_A_VARIANCE: variance of the lower truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_b_cdf(double x, double mu, double sigma, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_B_CDF evaluates the upper truncated Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the CDF.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+
|
|
|
+ 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;
|
|
|
+
|
|
|
+ return cdf;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void truncated_normal_b_cdf_values(int &n_data, double &mu, double &sigma,
|
|
|
+ double &b, double &x, double &fx)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 (-oo,B].
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, double &MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Output, double &SIGMA, the standard deviation of the distribution.
|
|
|
+//
|
|
|
+// Output, double &B, the upper truncation limit.
|
|
|
+//
|
|
|
+// Output, double &X, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_b_cdf_inv(double cdf, double mu, double sigma, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_B_CDF_INV inverts the upper truncated Normal CDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double CDF, the value of the CDF.
|
|
|
+// 0.0 <= CDF <= 1.0.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_B_CDF_INV - Fatal error!\n";
|
|
|
+ cerr << " 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_b_mean(double mu, double sigma, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_B_MEAN returns the mean of the upper truncated Normal
|
|
|
+// PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_b_moment(int order, double mu, double sigma, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_B_MOMENT: moments of the upper truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 11 September 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Reference:
|
|
|
+//
|
|
|
+// Phoebus Dhrymes,
|
|
|
+// Moments of Truncated Normal Distributions,
|
|
|
+// May 2005.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, int ORDER, the order of the moment.
|
|
|
+// 0 <= ORDER.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n";
|
|
|
+ cerr << " 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) {
|
|
|
+ cerr << "\n";
|
|
|
+ cerr << "TRUNCATED_NORMAL_B_MOMENT - Fatal error!\n";
|
|
|
+ cerr << " CDF((B-MU)/SIGMA) = 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_b_pdf(double x, double mu, double sigma, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_B_PDF evaluates the upper truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double X, the argument of the PDF.
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+
|
|
|
+ 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;
|
|
|
+
|
|
|
+ return pdf;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+void truncated_normal_b_pdf_values(int &n_data, double &mu, double &sigma,
|
|
|
+ double &b, double &x, double &fx)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// 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 GNU LGPL 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.
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input/output, int &N_DATA. The user sets N_DATA to 0 before the
|
|
|
+// first call. 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.
|
|
|
+//
|
|
|
+// Output, double &MU, the mean of the distribution.
|
|
|
+//
|
|
|
+// Output, double &SIGMA, the standard deviation of the distribution.
|
|
|
+//
|
|
|
+// Output, double &B, the upper truncation limit.
|
|
|
+//
|
|
|
+// Output, double &X, the argument of the function.
|
|
|
+//
|
|
|
+// Output, 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
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_b_sample(double mu, double sigma, double b, int &seed)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_B_SAMPLE samples the upper truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 21 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, double B, the upper truncation limit.
|
|
|
+//
|
|
|
+// Input/output, int &SEED, a seed for the random number
|
|
|
+// generator.
|
|
|
+//
|
|
|
+// Output, 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;
|
|
|
+}
|
|
|
+//****************************************************************************80
|
|
|
+
|
|
|
+double truncated_normal_b_variance(double mu, double sigma, double b)
|
|
|
+
|
|
|
+//****************************************************************************80
|
|
|
+//
|
|
|
+// Purpose:
|
|
|
+//
|
|
|
+// TRUNCATED_NORMAL_B_VARIANCE: variance of the upper truncated Normal PDF.
|
|
|
+//
|
|
|
+// Licensing:
|
|
|
+//
|
|
|
+// This code is distributed under the GNU LGPL license.
|
|
|
+//
|
|
|
+// Modified:
|
|
|
+//
|
|
|
+// 14 August 2013
|
|
|
+//
|
|
|
+// Author:
|
|
|
+//
|
|
|
+// John Burkardt
|
|
|
+//
|
|
|
+// Parameters:
|
|
|
+//
|
|
|
+// Input, double MU, SIGMA, the mean and standard deviation of the
|
|
|
+// parent Normal distribution.
|
|
|
+//
|
|
|
+// Input, 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;
|
|
|
+}
|
