tryrand3.cpp

随机数发生器C++写的

#define WANT_STREAM
#define WANT_MATH
#define WANT_TIME

#include "include.h"
#include "newran.h"
#include "tryrand.h"

#ifdef use_namespace
using namespace NEWRAN;
#endif

void SortAscending(Real* data, int max);
Real KS(Real* data, int n);
Real NormalDF(Real x);
double invchi95(int N);
double invchi99(int N);
void ChiSquaredTest(int* Observed, Real* Prob, int N, int n);
void TestBinomial(int N, Real p, int n);
void TestPoisson(Real mu, int n);
void TestNegativeBinomial(Real NX, Real p, int n);
void TestDiscreteGen(int N, Real* prob, int n);

inline Real square(Real x) { return x*x; }
inline Real cube(Real x) { return x*x*x; }


void test3(int n)
{
   cout << endl;

   // Do chi-squared tests to discrete data
   cout << "ChiSquared tests for discrete data" << endl;
   cout << "chisq should be less than 95% point in most cases" << endl;
   cout << "   and 99% point in almost all cases" << endl << endl;
   {
      Real p[] = { 0.05, 0.10, 0.05, 0.5, 0.01, 0.01, 0.03, 0.20, 0.05 };
      TestDiscreteGen(9, p, n);
   }

   {
      Real p[] = { 0.4, 0.2, 0.1, 0.05, 0.025, 0.0125, 0.00625, 0.00625, 0.2 };
      TestDiscreteGen(9, p, n);
   }


   TestNegativeBinomial(200.3, 0.05, n);
   TestNegativeBinomial(150.3, 0.15, n);
   TestNegativeBinomial(100.8, 0.18, n);
   TestNegativeBinomial(100.8, 1.22, n);
   TestNegativeBinomial(100.8, 9.0, n);
   TestNegativeBinomial(10.5, 0.18, n);
   TestNegativeBinomial(10.5, 1.22, n);
   TestNegativeBinomial(10.5, 9.0, n);
   TestNegativeBinomial(0.35, 0.18, n);
   TestNegativeBinomial(0.35, 1.22, n);
   TestNegativeBinomial(0.35, 9.0, n);

   TestBinomial(100, 0.45, n);
   TestBinomial(100, 0.25, n);
   TestBinomial(100, 0.02, n);
   TestBinomial(100, 0.01, n);
   TestBinomial(49, 0.60, n);
   TestBinomial(21, 0.70, n);
   TestBinomial(10, 0.90, n);
   TestBinomial(10, 0.25, n);
   TestBinomial(10, 0.10, n);

   TestPoisson(0.75, n);
   TestPoisson(4.3, n);
   TestPoisson(10, n);
   TestPoisson(100, n);

   Real* data = new Real[n];
   if (!data) Throw(Bad_alloc());

// Apply KS test to a variety of continuous distributions
//    - use cdf transform to convert to uniform

   cout << endl;
   cout << "Kolmogorov-Smirnoff tests for continuous distributions" << endl;
   cout << "25%, 5%, 1%, .1% upper points are 1.019, 1.358, 1.628, 1.950"
      << endl;
   cout << "5% lower point is 0.520" << endl;
   cout << "Values should be mostly less than 5% upper point" << endl;
   cout << "   and less than 1% point almost always" << endl << endl;

   {
      ChiSq X(1, 1.44);
      for (int i = 0; i < n; i++)
      {
         Real x = sqrt(X.Next());
         data[i] = NormalDF(x - 1.2) - NormalDF(-x - 1.2);
      }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      ChiSq X(4);
      for (int i = 0; i < n; i++)
         { Real x = 0.5 * X.Next(); data[i] = (1+x)*exp(-x); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      ChiSq X(2);
      for (int i = 0; i < n; i++) data[i] = exp(-0.5 * X.Next());
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Pareto X(0.5);
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = 1.0 / sqrt(x); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Pareto X(1.5);
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = 1.0 / (x * sqrt(x)); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Normal X;
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = NormalDF(x); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Normal N; SumRandom X = 10 + 5 * N;
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = NormalDF((x-10)/5); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Normal N; Cauchy C; MixedRandom X = N(0.9) + C(0.1);
      for (int i = 0; i < n; i++)
      {
         Real x = X.Next();
         data[i] = 0.9*NormalDF(x)+0.1*(atan(x)/3.141592654 + 0.5);
      }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Normal N; MixedRandom X = N(0.9) + (10*N)(0.1);
      for (int i = 0; i < n; i++)
      {
         Real x = X.Next();
         data[i] = 0.9*NormalDF(x)+0.1*NormalDF(x/10);
      }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Normal  X0; SumRandom X = X0 * 0.6 + X0 * 0.8;
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = NormalDF(x); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Normal X1;
      MixedRandom X = X1(0.2) + (X1 * 2.5 + 1.1)(0.35) + (X1 + 2.3)(0.45);
      for (int i = 0; i < n; i++)
      {
         Real x = X.Next();
         data[i] = 0.20 * NormalDF(x)
                 + 0.35 * NormalDF((x - 1.1) / 2.5)
                 + 0.45 * NormalDF(x - 2.3);
      }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Gamma X(0.5);
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = 2.0 * NormalDF(-sqrt(2 * x)); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Gamma X(3);
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = (1+x+0.5*x*x)*exp(-x); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Gamma X1(0.85); Gamma X2(2.15); SumRandom X = X1 + X2;
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = (1+x+0.5*x*x)*exp(-x); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Gamma X1(0.75); Gamma X2(0.25); SumRandom X = X1 + X2;
      for (int i = 0; i < n; i++) data[i] = exp(-X.Next());
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Gamma X(2);
      for (int i = 0; i < n; i++)
         { Real x = X.Next(); data[i] = (1+x)*exp(-x); }
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Exponential X;
      for (int i = 0; i < n; i++) data[i] = exp(-X.Next());
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Cauchy X;
      for (int i = 0; i < n; i++) data[i] = atan(X.Next())/3.141592654 + 0.5;
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Cauchy X0; SumRandom X = X0 * 0.3 + X0 * 0.7;
      for (int i = 0; i < n; i++) data[i] = atan(X.Next())/3.141592654 + 0.5;
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   {
      Uniform X;
      for (int i = 0; i < n; i++) data[i] = X.Next();
      cout << X.Name() << ":   "  << KS(data, n) << endl;
   }

   delete [] data;


}

/*************************** Kolmogorov Smirnov Test ************************/

// test the data in the array (length n) for being uniform (0,1)

Real KS(Real* data, int n)
{
   SortAscending(data, n);
   Real D = 0.0;
   for (int i = 0; i < n; i++)
   {
      Real d1 = (Real)(i+1) / (Real)n - data[i];
      Real d2 = data[i] - (Real)i / (Real)n;
      if (D < d1) D = d1; if (D < d2) D = d2;
   }
   return D * (sqrt((Real)n) + 0.12 + 0.11 / sqrt((Real)n));
}




/******************************** Quick sort ********************************/

// Quicksort.
// Essentially the method described in Sedgewick's algorithms in C++
// My version is still partially recursive, unlike Segewick's, but the
// smallest segment of each split is used in the recursion, so it should
// not overlead the stack.

// If the process does not seems to be converging an exception is thrown.


#define DoSimpleSort 17            // when to switch to insert sort
#define MaxDepth 50                // maximum recursion depth