home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 41 / IOPROG_41.ISO / soft / c++ / NUMCPP11.ZIP / rankdist.cpp < prev    next >
Encoding:
C/C++ Source or Header  |  1999-05-09  |  10.6 KB  |  476 lines
  1. //===================================================================
  2. // rankdist.cpp
  3. //
  4. // Version 1.1
  5. //
  6. // Written by:
  7. //   Brent Worden
  8. //   WordenWare
  9. //   email:  Brent@Worden.org
  10. //
  11. // Copyright (c) 1998-1999 WordenWare
  12. //
  13. // Created:  August 28, 1998
  14. // Revised:  April 10, 1999
  15. //===================================================================
  16.  
  17. #include <cmath>
  18.  
  19. #include "normdist.h"
  20. #include "rankdist.h"
  21. #include "vector.hpp"
  22.  
  23. NUM_BEGIN
  24.  
  25. // Function needed for Ansary-Bradley
  26. int start1(int n, double *f)
  27. {
  28.     int i;
  29.     int ret = n/2 + 1;
  30.     
  31.     for (i = 0; i < ret; ++i) {
  32.         f[i] = 2.0;
  33.     }
  34.     if(n % 2 == 0) {
  35.         f[ret-1] = 1.0;
  36.     }
  37.     
  38.     return ret;
  39. }
  40.  
  41. // Function needed for Ansary-Bradley
  42. int start2(int n, double *f)
  43. {
  44.     double a, b;
  45.     int i, j, nu, lt1, ndo;
  46.     int ret;
  47.     
  48.     nu = n - n % 2;
  49.     ndo = (lt1 = (ret = nu + 1) + 1) / 2;
  50.     a = 1.0;
  51.     b = 3.0;
  52.     for(i = 0, j = nu; i < ndo; ++i, --j) {
  53.         f[i] = a;
  54.         f[j] = a;
  55.         a += b;
  56.         b = 4.0 - b;
  57.     }
  58.     
  59.     if(nu != n){
  60.         nu = ndo + 1;
  61.         for (i = nu; i < ret; ++i) {
  62.             f[i] += 2.0;
  63.         }
  64.         f[lt1-1] = 2.0;
  65.         ret = lt1;
  66.     }
  67.     
  68.     return ret;
  69. }
  70.  
  71. // Function needed for Ansary-Bradley
  72. int frqadd(double *f1, int l1in, double *f2, int l2, int *nstart)
  73. {
  74.     int i1, i2, nxt, ret;
  75.     
  76.     for (i1 = *nstart-1, i2 = 0; i1 < l1in; ++i1, ++i2) {
  77.         f1[i1] += 2.0 * f2[i2];
  78.     }
  79.     nxt = l1in + 1;
  80.     ret = l2 + *nstart - 1;
  81.     for(i1 = nxt-1; i1 < ret; ++i1, ++i2){
  82.         f1[i1] = 2.0 * f2[i2];
  83.     }
  84.     ++(*nstart);
  85.     
  86.     return ret;
  87. }
  88.  
  89. int imply(double *f1, int l1in, int l1out, double *f2, int noff)
  90. {
  91.     double diff, sum;
  92.     int j2min, i2, j1, j2, i1, ndo, ret;
  93.     
  94.     i2 = 1 - noff - 1;
  95.     j1 = l1out;
  96.     j2 = l1out - noff;
  97.     ret = j2;
  98.     j2min = (j2 + 1) / 2;
  99.     ndo = (l1out + 1) / 2;
  100.     
  101.     for(i1 = 0; i1 < ndo; ++i1, ++i2, --j1){
  102.         if(i2 >= 0){
  103.             sum = f1[i1] + f2[i2];
  104.             f1[i1] = sum;
  105.         } else {
  106.             sum = f1[i1];
  107.         }
  108.         if(j2 >= j2min){
  109.             if(j1 > l1in){
  110.                 diff = sum;
  111.             } else {
  112.                 diff = sum - f1[j1-1];
  113.             }
  114.             f2[i1] = diff;
  115.             f2[j2-1] = diff;
  116.             --j2;
  117.         }
  118.         f1[j1-1] = sum;
  119.     }
  120.     
  121.     return ret;
  122. }
  123.  
  124. // Function needed for Ansary-Bradley
  125. void abmass(int m, int n, double *a1, int l1)
  126. {
  127.     int i1, i2, lres, mnow, symm, l1out, l2out, i, j, mm, nn;
  128.     int nc, mm1, ln1, nm1, nm2, ln3, ln2, n2b1, n2b2, ndo, astart;
  129.     double ai, *a2, *a3;
  130.     
  131.     a2 = new double[l1+1];
  132.     a3 = new double[l1+1];
  133.     
  134.     mm = ((m < n) ? m : n);
  135.     if (mm < 0) {
  136.         delete [] a2;
  137.         delete [] a3;
  138.         return;
  139.     }
  140.     i1 = (m + 1) / 2;
  141.     i2 = m / 2 + 1;
  142.     astart = i1 * i2;
  143.     nn = ((m > n) ? m : n);
  144.     
  145.     lres = mm * nn / 2 + 1;
  146.     if(l1 < lres){
  147.         delete [] a2;
  148.         delete [] a3;
  149.         return;
  150.     }
  151.     symm = (mm + nn) % 2 == 0;
  152.     
  153.     mm1 = mm - 1;
  154.     if(mm > 2){
  155.         nm1 = nn - 1;
  156.         nm2 = nn - 2;
  157.         mnow = 3;
  158.         nc = 3;
  159.         if(nn % 2 == 1){
  160.             n2b1 = 2;
  161.             n2b2 = 3;
  162.             ln1 = start1(nn, a1);
  163.             ln2 = start2(nm1, a2);
  164.             do {
  165.                 l1out = frqadd(a1, ln1, a2, ln2, &n2b1);
  166.                 ln1 += nn;
  167.                 ln3 = imply(a1, l1out, ln1, a3, nc);
  168.                 ++nc;
  169.                 if(mnow != mm){
  170.                     ++mnow;
  171.                     l2out = frqadd(a2, ln2, a3, ln3, &n2b2);
  172.                     ln2 += nm1;
  173.                     j = imply(a2, l2out, ln2, a3, nc);
  174.                     ++nc;
  175.                     if(mnow != mm){
  176.                         ++mnow;
  177.                     }
  178.                 }
  179.             } while(mnow != mm);
  180.         } else {
  181.             n2b1 = 3;
  182.             n2b2 = 2;
  183.             ln1 = start2(nn, a1);
  184.             ln3 = start2(nm2, a3);
  185.             ln2 = start1(nm1, a2);
  186.             do {
  187.                 l2out = frqadd(a2, ln2, a3, ln3, &n2b2);
  188.                 ln2 += nm1;
  189.                 j = imply(a2, l2out, ln2, a3, nc);
  190.                 ++nc;
  191.                 if(mnow != mm){
  192.                     ++mnow;
  193.                     l1out = frqadd(a1, ln1, a2, ln2, &n2b1);
  194.                     ln1 += nn;
  195.                     ln3 = imply(a1, l1out, ln1, a3, nc);
  196.                     ++nc;
  197.                     if(mnow != mm){
  198.                         ++mnow;
  199.                     }
  200.                 }
  201.             } while(mnow != mm);
  202.         }
  203.         if(symm){
  204.             delete [] a2;
  205.             delete [] a3;
  206.             return;
  207.         }
  208.         
  209.         for(i = (mm+3)/2, j = 1; i <= lres; ++i, ++j) {
  210.             if(i > ln1) {
  211.                 a1[i-1] = a2[j-1];
  212.             } else {
  213.                 a1[i-1] += a2[j-1];
  214.             }
  215.         }
  216.         
  217.         if(n < m){
  218.             delete [] a2;
  219.             delete [] a3;
  220.             return;
  221.         }
  222.     } else {
  223.         if(mm1 < 0){
  224.             a1[1-1] = 1.0;
  225.             delete [] a2;
  226.             delete [] a3;
  227.             return;
  228.         } else if(mm1 == 0){
  229.             ln1 = start1(nn, a1);
  230.         } else {
  231.             ln1 = start2(nn, a1);
  232.         }
  233.         if(symm || n > m){
  234.             delete [] a2;
  235.             delete [] a3;
  236.             return;
  237.         }
  238.     }
  239.     ndo = lres / 2;
  240.     for(i = 1, j = lres; i <= ndo; ++i, --j) {
  241.         ai = a1[i-1];
  242.         a1[i-1] = a1[j-1];
  243.         a1[j-1] = ai;
  244.     }
  245.     delete [] a2;
  246.     delete [] a3;
  247. }
  248.  
  249. NUMERICS_EXPORT double ansbradp(double x, int m, int n)
  250. {
  251.     int astart;
  252.     int i, nrows;
  253.     double sum;
  254.     int test = m, other = n;
  255.     double *a1;
  256.     double ret = 0.0;
  257.     int l1 = (m*n)/2+1;
  258.     
  259.     a1 = new double[l1+1];
  260.     
  261.     astart = m/2 * (m/2 + 1);
  262.     if(m % 2 == 1){
  263.         astart += (m+1)/2;
  264.     }
  265.     
  266.     abmass(m, n, a1, l1);
  267.     
  268.     nrows = m * n / 2 + 1;
  269.     sum = 0.0;
  270.     for (i = 0; i < nrows; ++i) {
  271.         sum += a1[i];
  272.         a1[i] = sum;
  273.     }
  274.     for (i = 0; i <= x-astart; ++i) {
  275.         a1[i] /= sum;
  276.     }
  277.     
  278.     ret = a1[(int)(x-astart)];
  279.     
  280.     delete [] a1;
  281.     
  282.     return ret;
  283. }
  284.  
  285. // Function needed for Wilcoxon
  286. int wilcox1(int x, int n, int k)
  287. {
  288.     int c, h, m, w, v, j, y;
  289.     
  290.     if(k == 1){
  291.         if(x < n) return x;
  292.         return n;
  293.     }
  294.     c = 1;
  295.     h = k-1;
  296.     m = n-1;
  297.     w = 0;
  298.     for(j = 1; j <= h; j++) c = c*(m-j+1)/j;
  299.     v = h*(m-h)+h*(h+1)/2;
  300.     y=x-h-1;
  301.     while(y >= v){
  302.         w += c;
  303.         c = (c*(m-h))/m;
  304.         m--;
  305.         v-=h;
  306.         y=y-h-1;
  307.     }
  308.     do {
  309.         w += wilcox1(y, m, h);
  310.         m--;
  311.         y = y-h-1;
  312.     } while(2*y >= h*(h+1));
  313.     
  314.     return w;
  315. }
  316.  
  317. NUMERICS_EXPORT double wilcoxonp(int x, int n)
  318. {
  319.     int c, k, u, v, w, xx;
  320.     bool a;
  321.     double m=n*(n+1.0)/4.0, p, cv;
  322.     
  323.     if(n <= 25){
  324.         if(x < 0) return 0.0;
  325.         m = n*(n+1)/2;
  326.         if(x >= m) return 1.0;
  327.         if(2*x > m){
  328.             xx = (int)(m)-x-1;
  329.             a = false;
  330.         } else {
  331.             xx = x;
  332.             a = true;
  333.         }
  334.         c = 1; k = 1; u = 1; v = n; w = 1;
  335.         while(xx >= v){
  336.             c = (c*(n-k+1))/k;
  337.             w += c;
  338.             k++;
  339.             u += k;
  340.             v += n+1-k;
  341.         }
  342.         do {
  343.             w += wilcox1(xx, n, k);
  344.             k++;
  345.             u += k;
  346.         } while(xx >= u);
  347.         if(a) p = ldexp((double)w, -(int)(n));
  348.         else p = 1.0 - ldexp((double)w, -(int)(n));
  349.     } else {
  350.         cv = n*(n+1.0)*(2.0*n+1.0)/24.0;
  351.         p = normalp(x+.5, m, cv);
  352.     }
  353.     
  354.     return p;
  355. }
  356.  
  357. NUMERICS_EXPORT void wilcoxonv(double p, int n, int *w0, int *w1)
  358. {
  359.     double m, v, wd;
  360.     int w;
  361.     
  362.     m = (double)(n*(n+1))/4.0;
  363.     v = (double)(n*(n+1)*(2.0*n+1))/24.0;
  364.     wd = normalv(p, m, v);
  365.     if(p < 0.0 || p > 1.0 || n < 1){
  366.         *w0 = *w1 = -1;
  367.     } else if(n > 25){
  368.         if(wd < 1.0) wd = 1.0;
  369.         else if(wd > 2.0*m) wd = 2.0*m-.5;
  370.         *w0 = (int)floor(wd);
  371.         *w1 = (int)ceil(wd);
  372.     } else {
  373.         w = (int)wd;
  374.         while(wilcoxonp(w, n) > p) w--;
  375.         while(wilcoxonp(w, n) <= p) w++;
  376.         if(w <= 1){
  377.             *w0 = w;
  378.         } else *w0 = w-1;
  379.         *w1 = w;
  380.     }
  381. }
  382.  
  383. NUMERICS_EXPORT double wilmannwhitp(int w, int m, int n)
  384. {
  385.     int u = w - m * (m + 1) / 2;
  386.     int minmn = ((m < n) ? m : n);
  387.     int mn1 = m*n+1, i;
  388.     int maxmn = ((m > n) ? m : n), n1 = maxmn + 1;
  389.     double ret = 0.0;
  390.     
  391.     if(u < 0){
  392.         return ret;
  393.     }
  394.     
  395.     if(minmn < 1){
  396.         return ret;
  397.     }
  398.     
  399.     Vector<double> freq(mn1);
  400.     
  401.     for(i = 1; i <= n1; ++i){
  402.         freq[i-1] = 1.0;
  403.     }
  404.     
  405.     if(minmn != 1){
  406.         double sum;
  407.         int in = maxmn, l, k, j;
  408.         Vector<double> work((mn1+1)/2+minmn);
  409.         
  410.         n1 = n1 + 1;
  411.         for(i = n1; i <= mn1; ++i){
  412.             freq[i-1] = 0.0;
  413.         }
  414.         
  415.         work[1-1] = 0.0;
  416.         for(i = 2; i <= minmn; ++i){
  417.             work[i-1] = 0.0;
  418.             in += maxmn;
  419.             n1 = in + 2;
  420.             l = 1 + in/2;
  421.             k = i;
  422.             for(j = 1; j <= l; ++j){
  423.                 ++k;
  424.                 --n1;
  425.                 sum = freq[j-1] + work[j-1];
  426.                 freq[j-1] = sum;
  427.                 work[k-1] = sum - freq[n1-1];
  428.                 freq[n1-1] = sum;
  429.             }
  430.         }
  431.         sum = 0.0;
  432.         for(i = 1; i <= mn1; ++i){
  433.             sum += freq[i-1];
  434.             freq[i-1] = sum;
  435.         }
  436.         for(i = 1; i <= u+1; ++i){
  437.             freq[i-1] /= sum;
  438.         }
  439.         ret = freq[u];
  440.     }
  441.     
  442.     return ret;
  443. }
  444.  
  445. NUMERICS_EXPORT void wilmannwhitv(double p, int m, int n, int *w0, int *w1)
  446. {
  447.     double mn, v, wd;
  448.     int w;
  449.     
  450.     if(p < 0.0 || p > 1.0 || n < 1){
  451.         *w0 = *w1 = -1;
  452.     } else {
  453.         mn = (double)(m*(m+n+1))/2.0;
  454.         v = (double)(m*n*(m+n+1))/12.0;
  455.         wd = normalv(p, mn, v);
  456.         w = (int)wd;
  457.         while(wilmannwhitp(w, m, n) > p) w--;
  458.         while(wilmannwhitp(w, m, n) <= p) w++;
  459.         if(w <= m*(m+1)/2){
  460.             *w0 = m*(m+1)/2;
  461.         } else {
  462.             *w0 = w-1;
  463.         }
  464.         *w1 = w;
  465.     }
  466. }
  467.  
  468. NUM_END
  469.  
  470. //===================================================================
  471. // Revision History
  472. //
  473. // Version 1.0 - 08/28/1998 - New.
  474. // Version 1.1 - 04/10/1999 - Added Numerics namespace.
  475. //===================================================================
  476.