tryrand3.cpp

D-ITG2.4源代码

static Real SortThreeDescending(Real* a, Real* b, Real* c);
static void MyQuickSortAscending(Real* first, Real* last, int depth);
static void InsertionSortAscending(Real* first, const int length, int guard);



static Real SortThreeDescending(Real* a, Real* b, Real* c)
{
   // sort *a, *b, *c; return *b; optimise for already sorted
   if (*a >= *b)
   {
      if (*b >= *c) return *b;
      else if (*a >= *c) { Real x = *c; *c = *b; *b = x; return x; }
      else { Real x = *a; *a = *c; *c = *b; *b = x; return x; }
   }
   else if (*c >= *b) { Real x = *c; *c = *a; *a = x; return *b; }
   else if (*a >= *c) { Real x = *a; *a = *b; *b = x; return x; }
   else { Real x = *c; *c = *a; *a = *b; *b = x; return x; }
}



void SortAscending(Real* data, int max)
{
   if (max > DoSimpleSort) MyQuickSortAscending(data, data + max - 1, 0);
   InsertionSortAscending(data, max, DoSimpleSort);
}

static void InsertionSortAscending(Real* first, const int length,
   int guard)
// guard gives the length of the sequence to scan to find first
// element (eg guard = length)
{
   if (length <= 1) return;

   // scan for first element
   Real* f = first; Real v = *f; Real* h = f;
   if (guard > length) guard = length; int i = guard - 1;
   while (i--) if (v > *(++f)) { v = *f; h = f; }
   *h = *first; *first = v;

   // do the sort
   i = length - 1; f = first;
   while (i--)
   {
      Real* g = f++; h = f; v = *h;
      while (*g > v) *h-- = *g--;
      *h = v;
   }
}

static void MyQuickSortAscending(Real* first, Real* last, int depth)
{
   for (;;)
   {
      const int length = last - first + 1;
      if (length < DoSimpleSort) return;
      if (depth++ > MaxDepth)
         Throw(Exception("QuickSortAscending fails"));
      Real* centre = first + length/2;
      const Real test = SortThreeDescending(last, centre, first);
      Real* f = first; Real* l = last;
      for (;;)
      {
         while (*(++f) < test) {}
         while (*(--l) > test) {}
         if (l <= f) break;
         const Real temp = *f; *f = *l; *l = temp;
      }
      if (f > centre) { MyQuickSortAscending(l+1, last, depth); last = f-1; }
      else { MyQuickSortAscending(first, f-1, depth); first = l+1; }
   }
}

Real NormalDF(Real x)
{
   // from Abramowitz and Stegun - accuracy 7.5E-8
   // accuracy is absolute; not relative
   // eventually will need a better method
   // but good enough here
   Real t = 1.0 / (1.0 + 0.2316419 * fabs(x));
   t = ( 0.319381530
     + (-0.356563782
     + ( 1.781477937
     + (-1.821255978
     +   1.330274429 * t) * t) * t) * t) * t;
   t = 0.3989422804014326779399461 * exp(-0.5 * x * x) * t;
   return (x < 0) ? t : 1.0 - t;
}

void ChiSquaredTest(int* Observed, Real* Prob, int N, int n)
{
   // go for at least two expected observations per cell
   // work in from ends

   if (N <= 0) { cout << "no categories" << endl; return; }
   if (n <= 0) { cout << "no data" << endl; return; }

   int O1 = 0; Real E1 = 0.0; int O2 = 0; Real E2 = 0.0;
   Real CS = 0.0; int df = 0; int i = 0; int Ni = N; Real ToGo = n;
   for (;;)
   {
      O1 += Observed[i]; Real e1 = n * Prob[i]; E1 += e1; ToGo -= e1;
      if (E1 >= 2.0 && ToGo + E2 >= 2.0)
         { CS += square(O1 - E1) / E1; df += 1; O1 = 0; E1 = 0.0; }

      if (i == Ni) break;
      ++i;

      O2 += Observed[Ni]; Real e2 = n * Prob[Ni]; E2 += e2; ToGo -= e2;
      if (E2 >= 2.0 && ToGo + E1 >= 2.0)
         { CS += square(O2 - E2) / E2; df += 1; O2 = 0; E2 = 0.0; }

      if (i == Ni) break;
      --Ni;
   }

   E1 += E2; O1 += O2;
   if (E1 > 0.0) { CS += square(O1 - E1) / E1; df += 1; }
   if (fabs(ToGo) >= 0.01) cout << "chi-squared program fails  - ";

   cout << "chisq = " << CS << "; df = " << (df-1)
      << "; 95% pt. = " << invchi95(df-1)
      << "; 99% pt. = " << invchi99(df-1) << endl;
}


void TestBinomial(int N, Real p, int n)
{
   Binomial X(N, p);
   Real q = 1.0 - p; Real ln_p = log(p); Real ln_q = log(q);
   int* obs = new int [N+1]; if (!obs) Throw(Bad_alloc());
   Real* prob = new Real [N+1]; if (!prob) Throw(Bad_alloc());
   int i;
   for (i = 0; i <= N; i++)
   {
      obs[i] = 0;
      prob[i] = exp(ln_gamma(N+1) - ln_gamma(i+1) - ln_gamma(N-i+1)
         + i * ln_p + (N-i) * ln_q);
   }
   for (i = 0; i < n; i++)
   {
      int b = (int)X.Next();
      if (b < 0 || b > N) Throw(Logic_error("Binomial error"));
      obs[b]++;
   }
   cout << "Binomial: "; ChiSquaredTest(obs, prob, N, n);

   delete [] obs; delete [] prob;
}

void TestPoisson(Real mu, int n)
{
   Poisson X(mu);
   Real ln_mu = log(mu);
   int N = (int)(20 + mu + 10 * sqrt(mu));         // set upper bound
   if (N > n)
   {
      cout << "Poisson: range too large" << endl;
      return;
   }
   int* obs = new int [N+1]; if (!obs) Throw(Bad_alloc());
   Real* prob = new Real [N+1]; if (!prob) Throw(Bad_alloc());
   int i;
   for (i = 0; i <= N; i++)
      { obs[i] = 0; prob[i] = exp(i * ln_mu - mu - ln_gamma(i+1)); }
   for (i = 0; i < n; i++)
   {
      int b = (int)(X.Next());
      if (b < 0 || b > N) Throw(Logic_error("Poisson error"));
      obs[b]++;
   }
   cout << "Poisson: "; ChiSquaredTest(obs, prob, N, n);

   delete [] obs; delete [] prob;
}

void TestNegativeBinomial(Real NX, Real P, int n)
{
   NegativeBinomial X(NX, P);
   Real Q = 1.0 + P; Real p = 1.0 / Q; Real q = 1.0 - p;
   Real ln_p = log(p); Real ln_q = log(q);
   Real mean = NX * P; Real var = mean * Q;
   int N = (int)(20 + mean + 100 * sqrt(var));         // set upper bound
      // won't be good enough for large P
   if (N > n)
   {
      cout << "NegativeBinomial: range too large" << endl;
      return;
   }
   int* obs = new int [N+1]; if (!obs) Throw(Bad_alloc());
   Real* prob = new Real [N+1]; if (!prob) Throw(Bad_alloc());
   int i;
   for (i = 0; i <= N; i++)
   {
      obs[i] = 0;
      prob[i] = exp(ln_gamma(NX+i) - ln_gamma(i+1) - ln_gamma(NX)
         + NX * ln_p + i * ln_q);
   }
   for (i = 0; i < n; i++)
   {
      int b = (int)X.Next();
      if (b < 0 || b > N) Throw(Logic_error("NegativeBinomial error"));
      obs[b]++;
   }
   cout << "NegativeBinomial: "; ChiSquaredTest(obs, prob, N, n);

   delete [] obs; delete [] prob;
}

void TestDiscreteGen(int N, Real* prob, int n)
{
   DiscreteGen X(N, prob);
   int* obs = new int [N]; if (!obs) Throw(Bad_alloc());
   int i;
   for (i = 0; i < N; i++) obs[i] = 0;
   for (i = 0; i < n; i++)
   {
      int b = (int)X.Next();
      if (b < 0 || b >= N) Throw(Logic_error("DiscreteGen error"));
      obs[b]++;
   }
   cout << "DiscreteGen: "; ChiSquaredTest(obs, prob, N-1, n);

   delete [] obs;
}

// Calculate 95% point of chi-squared distribution

double invchi95(int N)
// upper 95% point of chi-squared distribution
{
   if (N < 0) Throw(Logic_error("Error in invchi95 arg"));
   if (N < 30)
   {
      double Q[] = { 0, 3.841459, 5.991465, 7.814728, 9.487729, 11.0705,
         12.59159, 14.06714, 15.50731, 16.91898, 18.30704, 19.67506,
         21.02601, 22.36199, 23.68475, 24.99576, 26.2962, 27.58709,
         28.86928, 30.14351, 31.4104, 32.6705, 33.9244, 35.1725,
         36.4151, 37.6525, 38.8852, 40.1133, 41.3372, 42.5569 };
      return Q[N];
   }
   else
   {
      double A = 1.0/(4.5 * N); double H = (-0.0002 * 60)/N;
      double Q = N * cube(1 - A + (1.645 - H) * sqrt(A));
      return Q;
   }
}

// Calculate 99% point of chi-squared distribution

double invchi99(int N)
// upper 99% point of chi-squared distribution
{
   if (N < 0) Throw(Logic_error("Error in invchi99 arg"));
   if (N < 30)
   {
      double Q[] = { 0, 6.63490, 9.21034, 11.3449, 13.2767, 15.0863,
         16.8119, 18.4753, 20.0902, 21.6660, 23.2093, 24.7250,
         26.2170, 27.6883, 29.1413, 30.5779, 31.9999, 33.4087,
         34.8053, 36.1908, 37.5662, 38.9321, 40.2894, 41.6384,
         42.9798, 44.3141, 45.6417, 46.9630, 48.2782, 49.5879 };
      return Q[N];
   }
   else
   {
      double A = 1.0/(4.5 * N); double H = (0.0008 * 60)/N;
      double Q = N * cube(1 - A + (2.326 - H) * sqrt(A));
      return Q;
   }
}