• Main Page
  • Related Pages
  • Modules
  • Namespaces
  • Classes
  • Files
  • File List
  • File Members

src/SampleTest.cpp

00001 /*
00002  *   This file is part of the Standard Portable Library (SPL).
00003  *
00004  *   SPL is free software: you can redistribute it and/or modify
00005  *   it under the terms of the GNU General Public License as published by
00006  *   the Free Software Foundation, either version 3 of the License, or
00007  *   (at your option) any later version.
00008  *
00009  *   SPL is distributed in the hope that it will be useful,
00010  *   but WITHOUT ANY WARRANTY; without even the implied warranty of
00011  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00012  *   GNU General Public License for more details.
00013  *
00014  *   You should have received a copy of the GNU General Public License
00015  *   along with SPL.  If not, see <http://www.gnu.org/licenses/>.
00016  */
00017 #include <spl/math/SampleTest.h>
00018 
00019 static void _Sort(int N, double *ar)
00020 {
00021         for (int x = 0; x < N; x++)
00022         {
00023                 for(int y = x; y < N; y++)
00024                 {
00025                         if( ar[x] > ar[y] )
00026                         {
00027                                 double d = ar[x];
00028                                 ar[x] = ar[y];
00029                                 ar[y] = d;
00030                         }
00031                 }
00032         }
00033 }
00034 
00035 struct _SPair
00036 {
00037         double val;
00038         bool one;
00039 };
00040 
00041 static void _ArrayPairInsert(struct _SPair *a, const int alen, double d, bool one)
00042 {
00043 /*
00044 ' insert n value into an n dimensional array
00045 '
00046 ' called by:  RankSQV
00047 */
00048         for(int x = 0; x < alen; x++)
00049         {
00050                 if (d >= a[x].val)
00051                 {
00052                         //'
00053                         //' insert the value
00054                         //'
00055                         for (int y = alen-1; y > x; y--)
00056                         {
00057                                 a[y].val = a[y - 1].val;
00058                                 a[y].one = a[y - 1].one;
00059                         }
00060                         a[x].val = d;
00061                         a[x].one = one;
00062                         return;
00063                 }
00064         }
00065         throw new Exception("shouldn't get here");
00066 }
00067 
00068 double _AvgRank(struct _SPair *pairs, const int pairslen, int x)
00069 {
00070         double target = pairs[x].val;
00071 
00072         for(int i = 0; i < pairslen; i++)
00073         {
00074                 if (pairs[i].val == target)
00075                 {
00076                         // average the ranks of the duplicates
00077                         int N = 1;
00078                         double sum = i + 1;
00079                         while( i < pairslen )
00080                         {
00081                                 if( pairs[i].val != pairs[i + 1].val )
00082                                 {
00083                                         break;
00084                                 }
00085                                 N = N + 1;
00086                                 i = i + 1;
00087                                 sum = sum + i + 1;
00088                         }
00089                         return sum / N;
00090                 }
00091         }
00092 
00093         throw new Exception("Shouldn't get here");
00094 }
00095 
00096 static double _SortSumR1(Sample *a, Sample *b)
00097 {
00098         int x;
00099 
00100         double r = 0;
00101 
00102         int N = a->N() + b->N();
00103         struct _SPair *pairs = new _SPair[N];
00104         if ( NULL == pairs )
00105         {
00106                 throw OutOfMemoryException();
00107         }
00108 
00109         for (x = 0; x < N; x++)
00110         {
00111                 pairs[x].val = -1;
00112         }
00113 
00114         for (x = 0; x < a->N(); x++)
00115         {
00116                 _ArrayPairInsert( pairs, N, a->Item(x), true );
00117         }
00118 
00119         for (x = 0; x < b->N(); x++)
00120         {
00121                 _ArrayPairInsert( pairs, N, b->Item(x), false );
00122         }
00123 
00124         for (x = 0; x < N; x++)
00125         {
00126                 if (pairs[x].one)
00127                 {
00128                         r = r + _AvgRank(pairs, N, x);
00129                 }
00130         }
00131 
00132         delete[] pairs;
00133 
00134         return r;
00135 }
00136 
00137 static double _Rankits(int N, int rank)
00138 {
00139         static double table[][20] = {
00140         //              1               2               3               4               5               6               7               8               9                                       10              11              12              13              14              15              16              17              18              19              20
00141         {/*1*/  -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,                           -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00142         {/*2*/  -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,                           -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00143         {/*3*/  -0.846, 0,              0.846,  -777,   -777,   -777,   -777,   -777,   -777,                           -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00144         {/*4*/  -1.029, -0.297, 0.297,  1.029,  -777,   -777,   -777,   -777,   -777,                           -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00145         {/*5*/  -1.163, -0.495, 0,              0.495,  1.163,  -777,   -777,   -777,   -777,                           -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00146         {/*6*/  -1.267, -0.642, -0.202, 0.202,  0.642,  1.267,  -777,   -777,   -777,                           -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00147         {/*7*/  -1.352, -0.757, -0.353, 0,              0.353,  0.757,  1.352,  -777,   -777,           /*7*/   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00148         {/*8*/  -1.424, -0.852, -0.473, -0.153, 0.153,  0.473,  0.852,  1.424,  -777,           /*8*/   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00149         {/*9*/  -1.485, -0.932, -0.572, -0.275, 0,              0.275,  0.572,  0.932,  1.485,          /*9*/   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00150         {/*10*/ -1.539, -1.001, -0.656, -0.376, -0.123, 0.123,  0.376,  0.656,  1.001,          /*10*/  1.539,  -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00151         {/*11*/ -1.586, -1.062, -0.729, -0.462, -0.225, 0,              0.225,  0.462,  0.729,          /*11*/  1.062,  1.586,  -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00152         {/*12*/ -1.629, -1.116, -0.793, -0.537, -0.312, -0.103, 0.103,  0.312,  0.537,          /*12*/  0.793,  1.116,  1.629,  -777,   -777,   -777,   -777,   -777,   -777,   -777,   -777},
00153         {/*13*/ -1.668, -1.164, -0.85,  -0.603, -0.388, -0.19,  0,              0.19,   0.388,          /*13*/  0.603,  0.85,   1.164,  1.668,  -777,   -777,   -777,   -777,   -777,   -777,   -777},
00154         {/*14*/ -1.703, -1.208, -0.901, -0.662, -0.456, -0.267, -0.088, 0.088,  0.267,          /*14*/  0.456,  0.662,  0.901,  1.208,  1.703,  -777,   -777,   -777,   -777,   -777,   -777},
00155         {/*15*/ -1.736, -1.248, -0.948, -0.715, -0.516, -0.335, -0.165, 0,              0.165,          /*15*/  0.335,  0.516,  0.715,  0.948,  1.248,  1.736,  -777,   -777,   -777,   -777,   -777},
00156         {/*16*/ -1.766, -1.285, -0.99,  -0.763, -0.57,  -0.396, -0.234, -0.077, 0.077,          /*16*/  0.234,  0.396,  0.57,   0.763,  0.99,   1.285,  1.766,  -777,   -777,   -777,   -777},
00157         {/*17*/ -1.794, -1.319, -1.029, -0.807, -0.619, -0.451, -0.295, -0.146, 0,                      /*17*/  0.146,  0.295,  0.451,  0.619,  0.807,  1.029,  1.319,  1.794,  -777,   -777,   -777},
00158         {/*18*/ -1.82,  -1.35,  -1.066, -0.848, -0.665, -0.502, -0.351, -0.208, -0.069,         /*18*/  0.069,  0.208,  0.351,  0.502,  0.665,  0.848,  1.066,  1.35,   1.82,   -777,   -777},
00159         {/*19*/ -1.844, -1.38,  -1.099, -0.886, -0.707, -0.548, -0.402, -0.264, -0.131,         /*19*/  0,              -0.131, 0.264,  0.402,  0.548,  0.707,  0.886,  1.099,  1.38,   1.844,  -777},
00160         {/*20*/ -1.867, -1.408, -1.131, -0.921, -0.745, -0.59,  -0.448, -0.315, -0.187,         /*20*/  -0.062, 0.062,  0.187,  0.315,  0.448,  0.59,   0.745,  0.921,  1.131,  1.408,  1.867}
00161         //              1               2               3               4               5               6               7               8               9                                       10              11              12              13              14              15              16              17              18              19              20
00162         };
00163 //    Dim table(3 To 20, 1 To 20) As Double
00164 //
00165 //    'N , RANK1, RANK2, RANK3, RANK4, RANK5, RANK6, RANK7, RANK8, RANK9, RANK10, RANK11, RANK12, RANK13, RANK14, RANK15, RANK16, RANK17, RANK18, RANK19, RANK20
00166 //    '3,
00167 //    table(3, 1) = -0.846: table(3, 2) = 0#: table(3, 3) = 0.846
00168 //    table(4, 1) = -1.029: table(4, 2) = -0.297: table(4, 3) = 0.297: table(4, 4) = 1.029
00169 //    table(5, 1) = -1.163: table(5, 2) = -0.495: table(5, 3) = 0#: table(5, 4) = 0.495: table(5, 5) = 1.163
00170 //    table(6, 1) = -1.267: table(6, 2) = -0.642: table(6, 3) = -0.202: table(6, 4) = 0.202: table(6, 5) = 0.642: table(6, 6) = 1.267
00171 //    table(7, 1) = -1.352: table(7, 2) = -0.757: table(7, 3) = -0.353: table(7, 4) = 0#: table(7, 5) = 0.353: table(7, 6) = 0.757: table(7, 7) = 1.352
00172 //    table(8, 1) = -1.424: table(8, 2) = -0.852: table(8, 3) = -0.473: table(8, 4) = -0.153: table(8, 5) = 0.153: table(8, 6) = 0.473: table(8, 7) = 0.852: table(8, 8) = 1.424
00173 //    table(9, 1) = -1.485: table(9, 2) = -0.932: table(9, 3) = -0.572: table(9, 4) = -0.275: table(9, 5) = 0#: table(9, 6) = 0.275: table(9, 7) = 0.572: table(9, 8) = 0.932: table(9, 9) = 1.485
00174 //    table(10, 1) = -1.539: table(10, 2) = -1.001: table(10, 3) = -0.656: table(10, 4) = -0.376: table(10, 5) = -0.123: table(10, 6) = 0.123: table(10, 7) = 0.376: table(10, 8) = 0.656: table(10, 9) = 1.001: table(10, 10) = 1.539
00175 //    table(11, 1) = -1.586: table(11, 2) = -1.062: table(11, 3) = -0.729: table(11, 4) = -0.462: table(11, 5) = -0.225: table(11, 6) = 0#: table(11, 7) = 0.225: table(11, 8) = 0.462: table(11, 9) = 0.729: table(11, 10) = 1.062: table(11, 11) = 1.586
00176 //    table(12, 1) = -1.629: table(12, 2) = -1.116: table(12, 3) = -0.793: table(12, 4) = -0.537: table(12, 5) = -0.312: table(12, 6) = -0.103: table(12, 7) = 0.103: table(12, 8) = 0.312: table(12, 9) = 0.537: table(12, 10) = 0.793: table(12, 11) = 1.116: table(12, 12) = 1.629
00177 //    table(13, 1) = -1.668: table(13, 2) = -1.164: table(13, 3) = -0.85: table(13, 4) = -0.603: table(13, 5) = -0.388: table(13, 6) = -0.19: table(13, 7) = 0#: table(13, 8) = 0.19: table(13, 9) = 0.388: table(13, 10) = 0.603: table(13, 11) = 0.85: table(13, 12) = 1.164: table(13, 13) = 1.668
00178 //    table(14, 1) = -1.703: table(14, 2) = -1.208: table(14, 3) = -0.901: table(14, 4) = -0.662: table(14, 5) = -0.456: table(14, 6) = -0.267: table(14, 7) = -0.088: table(14, 8) = 0.088: table(14, 9) = 0.267: table(14, 10) = 0.456: table(14, 11) = 0.662: table(14, 12) = 0.901: table(14, 13) = 1.208: table(14, 14) = 1.703
00179 //    table(15, 1) = -1.736: table(15, 2) = -1.248: table(15, 3) = -0.948: table(15, 4) = -0.715: table(15, 5) = -0.516: table(15, 6) = -0.335: table(15, 7) = -0.165: table(15, 8) = 0#: table(15, 9) = 0.165: table(15, 10) = 0.335: table(15, 11) = 0.516: table(15, 12) = 0.715: table(15, 13) = 0.948: table(15, 14) = 1.248: table(15, 15) = 1.736
00180 //    table(16, 1) = -1.766: table(16, 2) = -1.285: table(16, 3) = -0.99: table(16, 4) = -0.763: table(16, 5) = -0.57: table(16, 6) = -0.396: table(16, 7) = -0.234: table(16, 8) = -0.077: table(16, 9) = 0.077: table(16, 10) = 0.234: table(16, 11) = 0.396: table(16, 12) = 0.57: table(16, 13) = 0.763: table(16, 14) = 0.99: table(16, 15) = 1.285: table(16, 16) = 1.766
00181 //    table(17, 1) = -1.794: table(17, 2) = -1.319: table(17, 3) = -1.029: table(17, 4) = -0.807: table(17, 5) = -0.619: table(17, 6) = -0.451: table(17, 7) = -0.295: table(17, 8) = -0.146: table(17, 9) = 0#: table(17, 10) = 0.146: table(17, 11) = 0.295: table(17, 12) = 0.451: table(17, 13) = 0.619: table(17, 14) = 0.807: table(17, 15) = 1.029: table(17, 16) = 1.319: table(17, 17) = 1.794
00182 //    table(18, 1) = -1.82: table(18, 2) = -1.35: table(18, 3) = -1.066: table(18, 4) = -0.848: table(18, 5) = -0.665: table(18, 6) = -0.502: table(18, 7) = -0.351: table(18, 8) = -0.208: table(18, 9) = -0.069: table(18, 10) = 0.069: table(18, 11) = 0.208: table(18, 12) = 0.351: table(18, 13) = 0.502: table(18, 14) = 0.665: table(18, 15) = 0.848: table(18, 16) = 1.066: table(18, 17) = 1.35: table(18, 18) = 1.82
00183 //    table(19, 1) = -1.844: table(19, 2) = -1.38: table(19, 3) = -1.099: table(19, 4) = -0.886: table(19, 5) = -0.707: table(19, 6) = -0.548: table(19, 7) = -0.402: table(19, 8) = -0.264: table(19, 9) = -0.131: table(19, 10) = 0#: table(19, 11) = 0.131: table(19, 12) = 0.264: table(19, 13) = 0.402: table(19, 14) = 0.548: table(19, 15) = 0.707: table(19, 16) = 0.886: table(19, 17) = 1.099: table(19, 18) = 1.38: table(19, 19) = 1.844
00184 //    '20,
00185 //    table(20, 1) = -1.867: table(20, 2) = -1.408: table(20, 3) = -1.131: table(20, 4) = -0.921: table(20, 5) = -0.745: table(20, 6) = -0.59: table(20, 7) = -0.448: table(20, 8) = -0.315: table(20, 9) = -0.187: table(20, 10) = -0.062: table(20, 11) = 0.062: table(20, 12) = 0.187: table(20, 13) = 0.315: table(20, 14) = 0.448: table(20, 15) = 0.59: table(20, 16) = 0.745: table(20, 17) = 0.921: table(20, 18) = 1.131: table(20, 19) = 1.408: table(20, 20) = 1.867
00186 //
00187         if ( N < 1 || N > 20 )
00188         {
00189                 throw new InvalidArgumentException("Rankit N must be in the range [1,20]");
00190         }
00191         if ( rank < 1 || rank > 20 )
00192         {
00193                 throw new InvalidArgumentException("Rankit rank  must be in the range [1,20]");
00194         }
00195         double prob = table[ rank-1 ][ N-1 ];
00196         ASSERT( prob != -777 && prob < 2 && prob > -2 );
00197         return prob;
00198 }
00199 
00200 static double _AvgRankits(struct _SPair *pairs, int N, int x)
00201 {
00202         ASSERT_MEM( pairs, N * sizeof(struct _SPair) );
00203         double sum = 0;
00204         int ncount = 0;
00205         ASSERT_PTR(&pairs[x]);
00206         double target = pairs[x].val;
00207 
00208         for (int i = 0; i < N; i++)
00209         {
00210                 if (pairs[i].val == target)
00211                 {
00212                         /* average the ranks of the duplicates*/
00213                         ncount = 1;
00214                         sum = _Rankits(N, i + 1);
00215                         while (i < N && pairs[i].val == pairs[i + 1].val)
00216                         {
00217                                 ncount = ncount + 1;
00218                                 i = i + 1;
00219                                 sum = sum + _Rankits(N, i + 1);
00220                         }
00221                         return sum / ncount;
00222                 }
00223         }
00224         throw new Exception("Internal error in _AvgRankits");
00225 }
00226 
00227 SampleTest::SampleTest(const Sample& s1, const Sample& s2)
00228 : m_s1(s1), m_s2(s2)
00229 {
00230 }
00231 
00232 SampleTest::SampleTest(List<double> *s1, List<double> *s2)
00233 : m_s1(s1), m_s2(s2)
00234 {
00235 }
00236 
00237 SampleTest::SampleTest(Vector<double> *s1, Vector<double> *s2)
00238 : m_s1(s1), m_s2(s2)
00239 {
00240 }
00241 
00242 SampleTest::SampleTest(const double *s1, const int s1count, const double *s2, const int s2count)
00243 : m_s1(s1, s1count), m_s2(s2, s2count)
00244 {
00245 }
00246 
00247 SampleTest::SampleTest(const SampleTest& t)
00248 : m_s1(t.m_s1), m_s2(t.m_s2)
00249 {
00250 }
00251 
00252 SampleTest::~SampleTest()
00253 {
00254 }
00255 
00256 SampleTest& SampleTest::operator =(const SampleTest& t)
00257 {
00258         m_s1 = t.m_s1;
00259         m_s2 = t.m_s2;
00260 
00261         return *this;
00262 }
00263 
00267 double SampleTest::StudentsT(double *prob)
00268 {
00269         double var1 = m_s1.AverageVariance();
00270         double var2 = m_s2.AverageVariance();
00271 
00272         int df = m_s1.N() + m_s2.N() - 2;
00273 
00274         double svar = ((m_s1.N() - 1) * var1 + (m_s2.N() - 1) * var2) / df;
00275 
00276         double t = (m_s1.Mean() - m_s2.Mean()) / Math::Sqrt(svar * (1.0 / m_s1.N() + 1.0 / m_s2.N()));
00277         if ( NULL != prob )
00278         {
00279                 *prob = Betai(0.5 * df, 0.5, df / (df + t * t));
00280         }
00281 
00282         return t;
00283 }
00284 
00288 double SampleTest::Betai(double a, double b, double x)
00289 {
00290     double bt;
00291     
00292     if( x < 0 || x > 1 )
00293         {
00294                 throw new InvalidArgumentException("x out of range in betai");
00295         }
00296     
00297         if( x == 0 || x == 1 )
00298         {
00299                 bt = 0;
00300         }
00301         else
00302         {
00303                 bt = Math::Exp(Gammaln(a + b) - Gammaln(a) - Gammaln(b) + a * Math::Log(x) + b * Math::Log(1 - x));
00304         }
00305 
00306         if( x < (a + 1) / (a + b + 2) )
00307         {
00308                 return bt * Betacf(a, b, x) / a;
00309         }
00310         else
00311         {
00312                 return 1 - bt * Betacf(b, a, 1 - x) / b;
00313         }
00314 }
00315 
00319 double SampleTest::Gammaln(double xx)
00320 {
00321         double cof[6];
00322 
00323         cof[0] = 76.1800917294715;
00324         cof[1] = -86.5053203294168;
00325         cof[2] = 24.0140982408309;
00326         cof[3] = -1.23173957245015;
00327         cof[4] = 1.20865097386618E-03;
00328         cof[5] = -5.395239384953E-06;
00329 
00330         double y = xx;
00331         double x = xx;
00332         double tmp = x + 5.5;
00333         tmp = tmp - (x + 0.5) * Math::Log(tmp);
00334         double ser = 1.00000000019001;
00335 
00336         for (int j = 0; j < 6; j++)
00337         {
00338                 y = y + 1;
00339                 ser = ser + cof[j] / y;
00340         }
00341         return -tmp + Math::Log(2.506628274631 * ser / x);
00342 }
00343 
00344 
00345 double SampleTest::Gamma(double d)
00346 {
00347         return Math::Exp(Gammaln(d));
00348 }
00349 
00353 double SampleTest::Betacf(double a, double b, double x)
00354 {
00355         const double FPMIN = 1E-30;
00356         const double EPS = 0.0000003;
00357         const int MAXIT = 200;
00358 
00359         double qab = a + b;
00360         double qap = a + 1;
00361         double qam = a - 1;
00362         double c = 1;
00363         double d = 1 - qab * x / qap;
00364         if( Math::Abs(d) < FPMIN )
00365         {
00366                 d = FPMIN;
00367         }
00368 
00369         d = 1 / d;
00370         double h = d;
00371         
00372         /* 100 is the maximum iterations */
00373         int m;
00374         for (m = 1; m <= MAXIT; m++ )
00375         {
00376                 int m2 = 2 * m;
00377                 double aa = m * (b - m) * x / ((qam + m2) * (a + m2));
00378                 d = 1 + aa * d;
00379                 if (Math::Abs(d) < FPMIN)
00380                 {
00381                         d = FPMIN;
00382                 }
00383                 c = 1 + aa / c;
00384                 if( Math::Abs(c) < FPMIN )
00385                 {
00386                         c = FPMIN;
00387                 }
00388                 d = 1 / d;
00389                 h = h * d * c;
00390                 aa = -(a + m) * (qab + m) * x / ((a + m2) * (qap + m2));
00391                 d = 1 + aa * d;
00392                 if (Math::Abs(d) < FPMIN)
00393                 {
00394                         d = FPMIN;
00395                 }
00396                 c = 1 + aa / c;
00397                 if (c < FPMIN)
00398                 {
00399                         c = FPMIN;
00400                 }
00401                 d = 1 / d;
00402                 double del = d * c;
00403                 h = h * del;
00404                 if (Math::Abs(del - 1) < EPS)
00405                 {
00406                         /* finished */
00407                         break;
00408                 }
00409         }
00410         if (m == MAXIT)
00411         {
00412                 throw new InvalidArgumentException("betacf a or b too big, or MAXIt too small");
00413         }
00414         return h;
00415 }
00416 
00421 double SampleTest::FTest(double *fstat)
00422 {
00423         double df1;
00424         double df2;
00425 
00426         double avgVar1 = 0;
00427         double avgVar2 = 0;
00428         double ep = 0;
00429         double s;
00430 
00431         int j;
00432 
00433         ASSERT(NULL != fstat);
00434 
00435         for( j = 0; j < m_s1.N(); j++ )
00436         {
00437                 s = m_s1.Item(j) - m_s1.Mean();
00438                 ep = ep + s;
00439                 avgVar1 = avgVar1 + s * s;
00440         }
00441         avgVar1 = (avgVar1 - ep * ep / m_s1.N()) / (m_s1.Mean() - 1);
00442 
00443         for( j = 0; j < m_s2.N(); j++ )
00444         {
00445                 s = m_s2.Item(j) - m_s2.Mean();
00446                 ep = ep + s;
00447                 avgVar2 = avgVar2 + s * s;
00448         }
00449         avgVar2 = (avgVar2 - ep * ep / m_s2.N()) / (m_s2.N() - 1);
00450 
00451         if (avgVar1 > avgVar2)
00452         {
00453                 if (avgVar2 == 0)
00454                 {
00455                         *fstat = 0;
00456                 }
00457                 else
00458                 {
00459                         *fstat = avgVar1 / avgVar2;
00460                 }
00461                 df1 = m_s1.N(); - 1;
00462                 df2 = m_s2.N() - 1;
00463         }
00464         else
00465         {
00466                 if (avgVar1 == 0)
00467                 {
00468                         *fstat = 0;
00469                 }
00470                 else
00471                 {
00472                         *fstat = avgVar2 / avgVar1;
00473                 }
00474                 df1 = m_s2.N(); - 1;
00475                 df2 = m_s1.N() - 1;
00476         }
00477 
00478         double prob = 2 * Betai(0.5 * df2, 0.5 * df1, df2 / (df2 + df1 * *fstat));
00479         if (prob > 1)
00480         {
00481                 prob = 2 - prob;
00482         }
00483         ASSERT(prob <= 1);
00484 
00485         return prob;
00486 }
00487 
00488 double SampleTest::ApproximateTTest(double *prob, double *df)
00489 {    
00490         ASSERT(NULL != prob && NULL != df);
00491 
00492         double var1 = m_s1.Variance();
00493         double var2 = m_s2.Variance();
00494 
00495         double t = (m_s1.Mean() - m_s2.Mean()) / Math::Sqrt(var1 / m_s1.N() + var2 / m_s2.N());
00496         *df = Math::Pow((var1 / m_s1.N() + var2 / m_s2.N()), 2) / (Math::Pow((var1 / m_s1.N()), 2) / (m_s1.N() - 1) + Math::Pow((var2 / m_s2.N()), 2) / (m_s2.N() - 1));
00497         *prob = Betai(0.5 * *df, 0.5 * *df, *df / (*df + t * t));
00498 
00499         return t;
00500 }
00501 
00502 double SampleTest::OneSampleTTest()
00503 {
00504     if( m_s1.Variance() == 0 )
00505         {
00507                 return 4;
00508         }
00509     else
00510         {
00511                 return (m_s1.Mean() - m_s2.Mean()) / (Math::Sqrt(m_s1.Variance() / m_s1.N()));
00512         }
00513 }
00514 
00515 //'{This program performs the Brown-Forsythe & Levene Tests for homogeneity of
00516 //'variances for the two group case.
00517 //'Levene and Brown-Forsythe tests for homogeneity of variances (HOV)
00518 //'A important assumption in the t-test for mean differences is that the variances
00519 //'in the two groups are equal (homogeneous). Two powerful and commonly used tests
00520 //'of this assumption are the Levene test and the Brown-Forsythe modification of
00521 //'this test. However, it is important to realize that (1) the homogeneity of
00522 //'variances assumption is usually not as crucial as other assumptions for the
00523 //'t-test for mean differences, in particular in the case of equal n, and (2) that
00524 //'the tests described below are not necessarily very robust themselves (e.g.,
00525 //'Glass and Hopkins, 1996, p. 436, call these tests "fatally flawed" see also the
00526 //'description of these tests below). If you are concerned about a violation of the
00527 //'HOV assumption, it is always advisable to repeat the key analyses using
00528 //'nonparametric methods.
00529 //'Levene 's test (homogeneity of variances): For each dependent variable, a t-test
00530 //'for mean differences is performed on the absolute deviations of values from the
00531 //'respective group means. If the Levene test is statistically significant, then
00532 //'the hypothesis of homogeneous variances should be rejected.
00533 //'Brown & Forsythe's test (homogeneity of variances): Recently, some authors
00534 //'(e.g., Glass and Hopkins, 1996) have called into question the power of the
00535 //'Levene test for unequal variances. Specifically, the absolute deviation (from
00536 //'the group means) scores can be expected to be highly skewed thus, the normality
00537 //'assumption for t-test for mean differences of those absolute deviation scores is
00538 //'usually violated. This poses a particular problem when there is unequal n in the
00539 //'two groups that are to be compared. A more robust test that is very similar to
00540 //'the Levene test has been proposed by Brown and Forsythe (1974). Instead of
00541 //'performing the t-test for mean differences on the deviations from the mean, one
00542 //'can perform the analysis on the deviations from the group medians. Olejnik and
00543 //'Algina (1987) have shown that this test will give quite accurate error rates
00544 //'even when the underlying distributions for the raw scores deviate significantly
00545 //'from the normal distribution. However, as Glass and Hopkins (1996, p. 436) have
00546 //'pointed out, both the Levene test as well as the Brown-Forsythe modification
00547 //'suffer from what those authors call a "fatal flaw," namely, that both tests
00548 //'themselves rely on the homogeneity of variances assumption (of the absolute
00549 //'deviations from the means or medians) and hence, it is not clear how robust
00550 //'these tests are themselves in the presence of significant variance heterogeneity
00551 //'and unequal n.
00552 //'Program written, modified, or edited at StatSoft, Inc.}
00553 
00556 double SampleTest::Levene(double *prob)
00557 {
00558         if( m_s1.N() == 0 || m_s2.N() == 0 )
00559         {
00560                 throw new Exception("No data in samples");
00561         }
00562         
00563         double *z1 = new double[m_s1.N()];
00564         double *z2 = new double[m_s2.N()];
00565         
00566         double med1 = m_s1.Median();
00567         double med2 = m_s2.Median();
00568 
00569         int j;
00570         double sum1 = 0;
00571         double zmean1 = 0;
00572         double zmean2 = 0;
00573 
00574         for( j = 0; j < m_s1.N(); j++)
00575         {
00576                 z1[j] = Math::Abs(m_s1.Item(j) - m_s1.Mean());
00577                 sum1 = sum1 + m_s1.Item(j);
00578         }
00579         zmean1 = sum1 / m_s1.N();
00580 
00581         sum1 = 0;
00582         for( j = 0; j < m_s2.N(); j++)
00583         {
00584                 z2[j] = Math::Abs(m_s2.Item(j));
00585                 sum1 = sum1 + m_s2.Item(j);
00586         }
00587         zmean2 = sum1 / m_s2.N();
00588         sum1 = 0;
00589 
00590         double sumsqr1 = 0;
00591         double sumsqr2 = 0;
00592         double sum2 = 0;
00593 
00594         for( j = 0; j < m_s1.N(); j++ )
00595         {
00596                 sum1 = sum1 + z1[j];
00597                 sumsqr1 = sumsqr1 + z1[j] * z1[j];
00598         }
00599 
00600         for( j = 0; j < m_s2.N(); j++)
00601         {
00602                 sum2 = sum2 + z2[j];
00603                 sumsqr2 = sumsqr2 + z2[j] * z2[j];
00604         }
00605 
00606         double sumsqr = 0;
00607         double SST = 0;
00608         double SSTRT = 0;
00609         double SSE = 0;
00610         double f = 0;
00611 
00612         sumsqr = (sum1 + sum2) * (sum1 + sum2);
00613         SST = sumsqr1 + sumsqr2 - sumsqr / (m_s1.N() + m_s2.N());
00614         SSTRT = (sum1 * sum1) / m_s1.N() + (sum2 * sum2) / m_s2.N() - sumsqr / (m_s1.N() + m_s2.N());
00615         SSE = SST - SSTRT;
00616         //degrees of freedom = Count + B.Count - 2
00617         f = SSTRT / (SSE / m_s1.N() + m_s2.N() - 2);
00618         *prob = 1 - FDistr(f, 1, m_s1.N() + m_s2.N() - 2);
00619 
00620         delete[] z1;
00621         delete[] z2;
00622         
00623         return f;
00624 }
00625 
00626 double SampleTest::LeveneAlternate(double *tprob)
00627 {
00628 /*
00629 '  This is the alternate x2 < x1 way.
00630 '
00631     ' 1. the absolute value of the difference between each replicate and the
00632     '    mean for its repective treatment is calculated.
00633     
00634     ' 2.  A 2-tailed t-test is performed between the test treatment absolute
00635     '     differences and the reference treatment absolute differences.
00636 */
00637         List<double> vals;    
00638     int i;
00639     
00640     for (i = 0; i < m_s1.N(); i++ )
00641         {
00642                 vals.Add( Math::Abs(m_s1.Item(i) - m_s1.Mean()) );
00643         }
00644     Sample ds1(&vals);
00645     
00646         vals.Clear();
00647     for (i = 0; i < m_s2.N(); i++ )
00648         {
00649                 vals.Add( Math::Abs(m_s2.Item(i) - m_s2.Mean()) );
00650         }
00651     Sample ds2(&vals);
00652     
00653         SampleTest tst(ds1, ds2);
00654 
00655         double t = tst.StudentsT(tprob);
00656         return Math::Abs(t);
00657 }
00658 
00659 double SampleTest::FDistr(double x, double v, double w)
00660 {
00661         return (Gamma((v + w) / 2)) / (Gamma(v / 2) * Gamma(w / 2)) * Math::Pow((v / w), (v / 2)) * Math::Pow(x, (v / 2)) * Math::Pow((1 + (v / w) * x), (-(v + w) / 2));
00662 }
00663 
00669 double SampleTest::WTestSmall()
00670 {
00671         const double *a = NULL;
00672         int N = m_s1.N() + m_s2.N();
00673         if (N > 50)
00674         {
00675                 throw new Exception("This version of the w-test only handles up to 50 data points");
00676         }
00677 
00678         static const double a3[] = {0.7071};
00679         static const double a4[] = { 0.6872, 0.1677, -1 };
00680         static const double a5[] = { 0.6646, 0.2413, -1 };
00681         static const double a6[] = { 0.6431, 0.2806, 0.0875, -1 };
00682         static const double a7[] = { 0.6233, 0.3031, 0.1401, -1 };
00683         static const double a8[] = { 0.6052, 0.3164, 0.1743, 0.0561, -1 };
00684         static const double a9[] = { 0.5888, 0.3244, 0.1976, 0.0947, -1 };;
00685         static const double a10[] = { 0.5739, 0.3291, 0.2141, 0.1224, 0.0399, -1 };
00686         static const double a11[] = { 0.5601, 0.3315, 0.226, 0.1429, 0.0695, -1 };
00687         static const double a12[] = { 0.5475, 0.3325, 0.2347, 0.1586, 0.0922, 0.0303, -1 };
00688         static const double a13[] = { 0.5359, 0.3325, 0.2412, 0.1707, 0.1099, 0.0539, -1 };
00689         static const double a14[] = { 0.5251,0.3318,0.246,0.1802,0.124,0.0727,0.024,-1 };
00690         static const double a15[] = { 0.515,0.3306,0.2495,0.1878,0.1353,0.088,0.0433,-1 };
00691         static const double a16[] = { 0.5056,0.329,0.2521,0.1939,0.1447,0.1005,0.0593,0.0196,-1 };
00692         static const double a17[] = { 0.4968,0.3273,0.254,0.1988,0.1524,0.1109,0.0725,0.0359,-1 };
00693         static const double a18[] = { 0.4886,0.3253,0.2553,0.2027,0.1587,0.1197,0.0837,0.0496,0.0163,-1 };
00694         static const double a19[] = { 0.4808,0.3232,0.2561,0.2059,0.1641,0.1271,0.0932,0.0612,0.0303,-1 };
00695         static const double a20[] = { 0.4734,0.3211,0.2565,0.2085,0.1686,0.1334,0.1013,0.0711,0.0422,0.014,-1 };
00696         static const double a21[] = { 0.4643,0.3185,0.2578,0.2119,0.1736,0.1399,0.1092,0.0804,0.053,0.0263,-1 };
00697         static const double a22[] = { 0.459,0.3156,0.2571,0.2131,0.1764,0.1443,0.115,0.0878,0.0618,0.0368,0.0122,-1 };
00698         static const double a23[] = { 0.4542,0.3126,0.2563,0.2139,0.1787,0.148,0.1201,0.0941,0.0696,0.0459,0.0228,-1 };
00699         static const double a24[] = { 0.4493,0.3098,0.2554,0.2145,0.1807,0.1512,0.1245,0.0997,0.0764,0.0539,0.0321,0.0107,-1 };
00700         static const double a25[] = { 0.445,0.3069,0.2543,0.2148,0.1822,0.1539,0.1283,0.1046,0.0823,0.061,0.0403,0.02,-1 };
00701         static const double a26[] = { 0.4407,0.3043,0.2533,0.2151,0.1836,0.1563,0.1316,0.1089,0.0876,0.0672,0.0476,0.0284,0.0094,-1 };
00702         static const double a27[] = { 0.4366,0.3018,0.2522,0.2152,0.1848,0.1584,0.1346,0.1128,0.0923,0.0728,0.054,0.0358,0.0178,-1 };
00703         static const double a28[] = { 0.4328,0.2992,0.251,0.2151,0.1857,0.1601,0.1372,0.1162,0.0965,0.0778,0.0598,0.0424,0.0253,0.0084,-1 };
00704         static const double a29[] = { 0.4291,0.2968,0.2499,0.215,0.1864,0.1616,0.1395,0.1192,0.1002,0.0822,0.065,0.0483,0.032,0.0159,-1 };
00705         static const double a30[] = { 0.4254,0.2944,0.2487,0.2148,0.187,0.163,0.1415,0.1219,
00706                 0.1036,0.0862,0.0697,0.0537,0.0381,0.0227,0.0076,-1 };
00707         static const double a31[] = { 0.422,0.2921,0.2475,0.2145,0.1874,0.1641,0.1433,0.1243,0.1066,0.0899,
00708                 0.0739,0.0585,0.0435,0.0289,0.0144,-1 };
00709         static const double a32[] = { 0.4188,0.2898,0.2463,0.2141,0.1878,0.1651,0.1449,0.1265,0.1093,0.0931,
00710                 0.0777,0.0629,0.0485,0.0344,0.0206,0.0068,-1 };
00711         static const double a33[] = { 0.4156,0.2876,0.2451,0.2137,0.188,0.166,0.1463,0.1284,0.1118,0.0961,0.0812,
00712                 0.0669,0.053,0.0395,0.0262,0.0131,-1 };
00713         static const double a34[] = { 0.4127,0.2854,0.2439,0.2132,0.1882,0.1667,0.1475,0.1301,0.114,0.0988,0.0844,
00714                 0.0706,0.0572,0.0441,0.0314,0.0187,0.0062,-1 };
00715         static const double a35[] = { 0.4096,0.2834,0.2427,0.2127,0.1883,0.1673,0.1487,0.1317,0.116,0.1013,0.0873,
00716                 0.0739,0.061,0.0484,0.0361,0.0239,0.0119,-1 };
00717         static const double a36[] = { 0.4068,0.2813,0.2415,0.2121,0.1883,0.1678,0.1496,0.1331,0.1179,0.1036,0.09,
00718                 0.077,0.0645,0.0523,0.0404,0.0287,0.0172,0.0057,-1 };
00719         static const double a37[] = { 0.404,0.2794,0.2403,0.2116,0.1883,0.1683,0.1505,0.1344,0.1196,0.1056,0.0924,
00720                 0.0798,0.0677,0.0559,0.0444,0.0331,0.022,0.011,-1 };
00721         static const double a38[] = { 0.4015,0.2774,0.2391,0.211,0.1881,0.1686,0.1513,0.1356,0.1211,0.1075,0.0947,
00722                 0.0824,0.0706,0.0592,0.0481,0.0372,0.0264,0.0158,0.0053,-1 };
00723         static const double a39[] = { 0.3989,0.2755,0.238,0.2104,0.188,0.1689,0.152,0.1366,0.1225,0.1092,0.0967,
00724                 0.0848,0.0733,0.0622,0.0515,0.0409,0.0305,0.0203,0.0101,-1 };
00725         static const double a40[] = { 0.3964,0.2737,0.2368,0.2098,0.1878,0.1691,0.1526,0.1376,0.1237,0.1108,
00726                 0.0986,0.087,0.0759,0.0651,0.0546,0.0444,0.0343,0.0244,0.0146,0.0049,-1 };
00727         static const double a41[] = { 0.394,0.2719,0.2357,0.2091,0.1876,0.1693,0.1531,0.1384,0.1249,0.1123,
00728                 0.1004,0.0891,0.0782,0.0677,0.0575,0.0476,0.0379,0.0283,0.0188,0.0094,-1 };
00729         static const double a42[] = { 0.3917,0.2701,0.2345,0.2085,0.1874,0.1694,0.1535,0.1392,0.1259,
00730                 0.1136,0.102,0.0909,0.0804,0.0701,0.0602,0.0506,0.0411,0.0318,0.0227,0.0136,0.0045,-1 };
00731         static const double a43[] = { 0.3894,0.2684,0.2334,0.2078,0.1871,0.1695,0.1539,0.1398,0.1269,0.1149,
00732                 0.1035,0.0927,0.0824,0.0724,0.0628,0.0534,0.0442,0.0352,0.0263,0.0175,0.0087,-1 };
00733         static const double a44[] = { 0.3872,0.2667,0.2323,0.2072,0.1868,0.1695,0.1542,0.1405,0.1278,0.116,
00734                 0.1049,0.0943,0.0842,0.0745,0.0651,0.056,0.0471,0.0383,0.0296,0.0211,0.0126,0.0042,-1 };
00735         static const double a45[] = { 0.385,0.2651,0.2313,0.2065,0.1865,0.1695,0.1545,0.141,0.1286,0.117,0.1062,
00736                 0.0959,0.086,0.0765,0.0673,0.0584,0.0497,0.0412,0.0328,0.0245,0.0163,0.0081,-1 };
00737         static const double a46[] = { 0.383,0.2635,0.2302,0.2058,0.1862,0.1695,0.1548,0.1415,0.1293,0.118,0.1073,
00738                 0.0972,0.0876,0.0783,0.0694,0.0607,0.0522,0.0439,0.0357,0.0277,0.0197,0.0118,0.0039,-1 };
00739         static const double a47[] = { 0.3808,0.262,0.2291,0.2052,0.1859,0.1695,0.155,0.142,0.13,0.1189,0.1085,0.0986,
00740                 0.0892,0.0801,0.0713,0.0628,0.0546,0.0465,0.0385,0.0307,0.0229,0.0153,0.0076,-1 };
00741         static const double a48[] = { 0.3789,0.2604,0.2281,0.2045,0.1855,0.1693,0.1551,0.1423,0.1306,0.1197,0.1095,
00742                 0.0998,0.0906,0.0817,0.0731,0.0648,0.0568,0.0489,0.0411,0.0335,0.0259,0.0185,0.0111,0.0037,-1 };
00743         static const double a49[] = { 0.377,0.2589,0.2271,0.2038,0.1851,0.1692,0.1553,0.1427,0.1312,0.1205,0.1105,0.101,
00744                 0.0919,0.0832,0.0748,0.0667,0.0588,0.0511,0.0436,0.0361,0.0288,0.0215,0.0143,0.0071,-1 };
00745         static const double a50[] = { 0.3751,0.2574,0.226,0.2032,0.1847,0.1691,0.1554,0.143,0.1317,0.1212,0.1113,
00746         0.102,0.0932,0.0846,0.0764,0.0685,0.0608,0.0532,0.0459,0.0386,0.0314,0.0244,0.0174,0.0104,0.0035,-1 };
00747 
00748         switch ( N )
00749         {
00750         case 3:
00751                 a = &a3[0];
00752                 break;
00753         case 4:
00754                 a = &a4[0];
00755                 break;
00756         case 5:
00757                 a = &a5[0];
00758                 break;
00759         case 6:
00760                 a = &a6[0];
00761                 break;
00762         case 7:
00763                 a = &a7[0];
00764                 break;
00765         case 8:
00766                 a = &a8[0];
00767                 break;
00768         case 9:
00769                 a = &a9[0];
00770                 break;
00771         case 10:
00772                 a = &a10[0];
00773                 break;
00774         case 11:
00775                 a = &a11[0];
00776                 break;
00777         case 12:
00778                 a = &a12[0];
00779                 break;
00780         case 13:
00781                 a = &a13[0];
00782                 break;
00783         case 14:
00784                 a = &a14[0];
00785                 break;
00786         case 15:
00787                 a = &a15[0];
00788                 break;
00789         case 16:
00790                 a = &a16[0];
00791                 break;
00792         case 17:
00793                 a = &a17[0];
00794                 break;
00795         case 18:
00796                 a = &a18[0];
00797                 break;
00798         case 19:
00799                 a = &a19[0];
00800                 break;
00801         case 20:
00802                 a = &a20[0];
00803                 break;
00804         case 21:
00805                 a = &a21[0];
00806                 break;
00807         case 22:
00808                 a = &a22[0];
00809                 break;
00810         case 23:
00811                 a = &a23[0];
00812                 break;
00813         case 24:
00814                 a = &a24[0];
00815                 break;
00816         case 25:
00817                 a = &a25[0];
00818                 break;
00819         case 26:
00820                 a = &a26[0];
00821                 break;
00822         case 27:
00823                 a = &a27[0];
00824                 break;
00825         case 28:
00826                 a = &a28[0];
00827                 break;
00828         case 29:
00829                 a = &a29[0];
00830                 break;
00831         case 30:
00832                 a = &a30[0];
00833                 break;
00834         case 31:
00835                 a = &a31[0];
00836                 break;
00837         case 32:
00838                 a = &a32[0];
00839                 break;
00840         case 33:
00841                 a = &a33[0];
00842                 break;
00843         case 34:
00844                 a = &a34[0];
00845                 break;
00846         case 35:
00847                 a = &a35[0];
00848                 break;
00849         case 36:
00850                 a = &a36[0];
00851                 break;
00852         case 37:
00853                 a = &a37[0];
00854                 break;
00855         case 38:
00856                 a = &a38[0];
00857                 break;
00858         case 39:
00859                 a = &a39[0];
00860                 break;
00861         case 40:
00862                 a = &a40[0];
00863                 break;
00864         case 41:
00865                 a = &a41[0];
00866                 break;
00867         case 42:
00868                 a = &a42[0];
00869                 break;
00870         case 43:
00871                 a = &a43[0];
00872                 break;
00873         case 44:
00874                 a = &a44[0];
00875                 break;
00876         case 45:
00877                 a = &a45[0];
00878                 break;
00879         case 46:
00880                 a = &a46[0];
00881                 break;
00882         case 47:
00883                 a = &a47[0];
00884                 break;
00885         case 48:
00886                 a = &a48[0];
00887                 break;
00888         case 49:
00889                 a = &a49[0];
00890                 break;
00891         case 50:
00892                 a = &a50[0];
00893                 break;
00894         default:
00895                 throw new Exception("Unsupported value of N in w-test small");
00896         }
00897 
00898         double sumb = 0;
00899         double sumx = 0;
00900         double sumx2 = 0;
00901 
00902         int i;
00903         int k = N / 2;
00904 
00905         double *x1 = new double[N];
00906 
00907         // put the data from both sets into an array
00908         for (i = 0; i < m_s1.N(); i++)
00909         {
00910                 x1[i] = m_s1.Item(i) - m_s1.Mean();
00911         }
00912 
00913         i = m_s1.N();
00914         for (int j = 0; j < m_s2.N(); j++)
00915         {
00916                 x1[i] = m_s2.Item(j) - m_s2.Mean();
00917                 i++;
00918         }
00919 
00920         _Sort(N, x1);
00921         ASSERT_MEM( x1, N * sizeof(double) );
00922 
00923         for( i = 0; i < N; i++ )
00924         {
00925                 sumx = sumx + x1[i];
00926         }
00927         sumx = sumx / N;
00928 
00929         for (i = 0; i < N; i++)
00930         {
00931                 sumx2 = sumx2 + (x1[i] - sumx) * (x1[i] - sumx);
00932         }
00933 
00934         if (((N - 1) % 2) == 0 )
00935         {
00936                 for(i = 0; i < k; i++)
00937                 {
00938                         sumb = sumb + a[i] * (x1[N - (i + 1)] - x1[i]);
00939                 }
00940         }
00941         else
00942         {
00943                 for(i = 0; i <= k; i++)
00944                 {
00945                         sumb = sumb + a[i] * (x1[N - (i + 1)] - x1[i]);
00946                 }
00947         }
00948 
00949         delete [] x1;
00950 
00951         if (sumx2 == 0)
00952         {
00953                 return 0;
00954         }
00955         else
00956         {
00957                 return sumb * sumb / sumx2;
00958         }
00959 }
00960 
00963 double SampleTest::UStat()
00964 {    
00965         int N = m_s1.N() + m_s2.N();       
00966         double R1 = _SortSumR1(&m_s1, &m_s2);
00967 
00968         return m_s1.N() * m_s2.N() + (m_s1.N() * (m_s1.N() + 1)) / 2 - R1;
00969 }
00970 
00973 double SampleTest::UStatAlternate()
00974 {
00975         int N = m_s1.N() + m_s2.N();       
00976         double R1 = _SortSumR1(&m_s2, &m_s1);
00977         return m_s1.N() * m_s2.N() + (m_s2.N() * (m_s2.N() + 1)) / 2 - R1;
00978 }
00979 
00980 double SampleTest::WTestCriticalValue(int N, double alpha)
00981 {
00982         int category;
00983 
00984         if ( alpha <= 0.01 )
00985         {
00986                 category = 1;
00987         }
00988         else if ( alpha <= 0.02 )
00989         {
00990                 category = 2;
00991         }
00992         else if ( alpha <= 0.05 )
00993         {
00994                 category = 3;
00995         }
00996         else if ( alpha <= 0.1 )
00997         {
00998                 category = 4;
00999         }
01000         else
01001         {
01002                 category = 5;
01003         }
01004     
01005         static double table[][50] = {
01006     //A0_01, A0_02, A0_05, A0_10, A0_50
01007         /*3*/{0.753,    0.756,  0.767,  0.789,  0.959},
01008     /*4*/{0.687,        0.707,  0.748,  0.792,  0.935},
01009     /*5*/{0.686,        0.715,  0.762,  0.806,  0.927},
01010     /*6*/{0.713,        0.743,  0.788,  0.826,  0.927},
01011     /*7*/{0.73,         0.76,   0.803,  0.838,  0.928},
01012     /*8*/{0.749,        0.778,  0.818,  0.851,  0.932},
01013     /*9*/{0.764,        0.791,  0.829,  0.859,  0.935},
01014     /*10*/{0.781,       0.806,  0.842,  0.869,  0.938},
01015     /*11*/{0.792,       0.817,  0.85,   0.876,  0.94},
01016     /*12*/{0.805,       0.828,  0.859,  0.883,  0.943},
01017     /*13*/{0.814,       0.837,  0.866,  0.889,  0.945},
01018     /*14*/{0.825,       0.846,  0.874,  0.895,  0.947},
01019     /*15*/{0.835,       0.855,  0.881,  0.901,  0.95},
01020     /*16*/{0.844,       0.863,  0.887,  0.906,  0.952},
01021     /*17*/{0.851,       0.869,  0.892,  0.91,   0.954},
01022     /*18*/{0.858,       0.874,  0.897,  0.914,  0.956},
01023     /*19*/{0.863,       0.879,  0.901,  0.917,  0.957},
01024     /*20*/{0.868,       0.884,  0.905,  0.92,   0.959},
01025     /*21*/{0.873,       0.888,  0.908,  0.923,  0.96},
01026               {0.878,   0.892,  0.911,  0.926,  0.961},
01027           {0.881,       0.895,  0.914,  0.928,  0.962},
01028           {0.884,       0.898,  0.916,  0.93,   0.963},
01029           {0.888,       0.901,  0.918,  0.931,  0.964},
01030         /*26*/
01031           {0.891,       0.904,  0.92,   0.933,  0.965},
01032           {0.894,       0.906,  0.923,  0.935,  0.965},
01033           {0.896,       0.908,  0.924,  0.936,  0.966},
01034           {0.898,       0.91,   0.926,  0.937,  0.966},
01035           {0.9,         0.912,  0.927,  0.939,  0.967},
01036         /*31*/
01037           {0.902,       0.914,  0.929,  0.94,   0.967},
01038           {0.904,       0.915,  0.93,   0.941,  0.968},
01039           {0.906,       0.917,  0.931,  0.942,  0.968},
01040           {0.908,       0.919,  0.933,  0.943,  0.969},
01041           {0.91,        0.92,   0.934,  0.944,  0.969},
01042         /*36*/
01043           {0.912,       0.922,  0.935,  0.945,  0.97},
01044           {0.914,       0.924,  0.936,  0.946,  0.97},
01045           {0.916,       0.925,  0.938,  0.947,  0.971},
01046           {0.917,       0.927,  0.939,  0.948,  0.971},
01047           {0.919,       0.928,  0.94,   0.949,  0.972},
01048         /*41*/
01049           {0.92,        0.929,  0.941,  0.95,   0.972},
01050           {0.922,       0.93,   0.942,  0.951,  0.972},
01051           {0.923,       0.932,  0.943,  0.951,  0.973},
01052           {0.924,       0.933,  0.944,  0.952,  0.973},
01053           {0.926,       0.934,  0.945,  0.953,  0.973},
01054         /*46*/
01055           {0.927,       0.935,  0.945,  0.953,  0.974},
01056           {0.928,       0.936,  0.946,  0.954,  0.974},
01057           {0.929,       0.937,  0.947,  0.954,  0.974},
01058           {0.929,       0.937,  0.947,  0.955,  0.974},
01059           {0.93,        0.938,  0.947,  0.955,  0.974}
01060         };
01061         
01062         return table[category][N-1];
01063 }
01064 
01065 double SampleTest::TCriticalValueAlternate(int N, double alpha)
01066 {
01067         int alphacat;
01068 
01069         if ( N < 1 || N > 50 )
01070         {
01071                 throw new InvalidArgumentException("N for critical value must in the range [1,50]");
01072         }
01073 
01074         // get a table column for the alpha level
01075         if (alpha <= 0.025)
01076         {
01077                 alphacat = 3;
01078         }
01079         else if( alpha <= 0.05 )
01080         {
01081                 alphacat = 2;
01082         }
01083         else if (alpha <= 0.1)
01084         {
01085                 alphacat = 1;
01086         }
01087         else
01088         {
01089                 alphacat = 0;
01090         }
01091         
01092         static double crtab[][50] = {
01093                 //A0_20         A0_10   A0_05   A_025
01094     /*1*/{1.376,        3,              6.31,   12.796},
01095     /*2*/{1.061,        2,              2.92,   4.303},
01096     /*3*/{0.978,        2,              2.353,  3.182},
01097     /*4*/{0.941,        2,              2.132,  2.776},
01098     /*5*/{0.92,         1,              2.015,  2.571},
01099     /*6*/{0.906,        1,              1.943,  2.447},
01100     /*7*/{0.896,        1,              1.895,  2.365},
01101     /*8*/{0.889,        1,              1.86,   2.306},
01102     /*9*/{0.883,        1,              1.833,  2.262},
01103     /*10*/{0.879,       1,              1.812,  2.228},
01104     /*11*/{0.876,       1,              1.796,  2.201},
01105     /*12*/{0.873,       1,              1.782,  2.179},
01106     /*13*/{0.87,        1,              1.771,  2.16},
01107     /*14*/{0.868,       1,              1.761,  2.145},
01108     /*15*/{0.866,       1,              1.753,  2.131},
01109     /*16*/{0.865,       1,              1.746,  2.12},
01110     /*17*/{0.863,       1,              1.74,   2.11},
01111     /*18*/{0.862,       1,              1.734,  2.101},
01112     /*19*/{0.861,       1,              1.729,  2.093},
01113     /*20*/{0.86,        1,              1.725,  2.086},
01114     /*21*/{0.859,       1,              1.721,  2.08},
01115     /*22*/{0.858,       1,              1.717,  2.074},
01116     /*23*/{0.858,       1,              1.714,  2.069},
01117     /*24*/{0.857,       1,              1.711,  2.064},
01118     /*25*/{0.856,       1,              1.708,  2.06},
01119     /*26*/{0.856,       1,              1.71,   2.06},
01120     /*27*/{0.855,       1,              1.7,    2.05},
01121     /*28*/{0.855,       1,              1.7,    2.05},
01122     /*29*/{0.854,       1,              1.7,    2.04},
01123     /*30*/{0.854,       1,              1.7,    2.04},
01124         
01125               {0.851,   1,              1.68,   2.02},
01126               {0.851,   1,              1.68,   2.02},
01127               {0.851,   1,              1.68,   2.02},
01128               {0.851,   1,              1.68,   2.02},
01129               {0.851,   1,              1.68,   2.02},
01130               {0.851,   1,              1.68,   2.02},
01131               {0.851,   1,              1.68,   2.02},
01132               {0.851,   1,              1.68,   2.02},
01133               {0.851,   1,              1.68,   2.02},
01134         /*40*/{0.851,   1,              1.68,   2.02},
01135                 
01136               {0.848,   1,              1.67,   2},
01137               {0.848,   1,              1.67,   2},
01138               {0.848,   1,              1.67,   2},
01139               {0.848,   1,              1.67,   2},
01140               {0.848,   1,              1.67,   2},
01141               {0.848,   1,              1.67,   2},
01142               {0.848,   1,              1.67,   2},
01143               {0.848,   1,              1.67,   2},
01144               {0.848,   1,              1.67,   2},
01145         /*50*/{0.848,   1,              1.67,   2}
01146         };      
01147     
01148     return crtab[alphacat][N-1];
01149 }
01150 
01153 double SampleTest::TCriticalValue(int N, double alpha)
01154 {
01155         int alphacat;
01156 
01157         if ( N < 1 || N > 50 )
01158         {
01159                 throw new InvalidArgumentException("N for critical value must in the range [1,50]");
01160         }
01161         
01162         // get a table column for the alpha level
01163         if (alpha <= 0.025)
01164         {
01165                 alphacat = 3;
01166         }
01167         else if( alpha <= 0.05 )
01168         {
01169                 alphacat = 2;
01170         }
01171         else if (alpha <= 0.1)
01172         {
01173                 alphacat = 1;
01174         }
01175         else
01176         {
01177                 alphacat = 0;
01178         }
01179         
01180         static double crtab[][50] = {
01181                 //A0_20         A0_10   A0_05   A_025
01182         /*1*/{1.376,    3.08,   6.31,   12.796},
01183     /*2*/{1.061,        1.89,   2.92,   4.303},
01184     /*3*/{0.978,        1.64,   2.353,  3.182},
01185     /*4*/{0.941,        1.53,   2.132,  2.776},
01186     /*5*/{0.92,         1.48,   2.015,  2.571},
01187     /*6*/{0.906,        1.44,   1.943,  2.447},
01188     /*7*/{0.896,        1.42,   1.895,  2.365},
01189     /*8*/{0.889,        1.4,    1.86,   2.306},
01190     /*9*/{0.883,        1.38,   1.833,  2.262},
01191     /*10*/{0.879,       1.37,   1.812,  2.228},
01192     /*11*/{0.8,         1.36,   1.796,  2.201},
01193     /*12*/{0.873,       1.36,   1.782,  2.179},
01194     /*13*/{0.87,        1.35,   1.771,  2.16},
01195     /*14*/{0.868,       1.34,   1.761,  2.145},
01196     /*15*/{0.866,       1.34,   1.753,  2.131},
01197     /*16*/{0.865,       1.34,   1.746,  2.12},
01198     /*17*/{0.863,       1.33,   1.74,   2.11},
01199     /*18*/{0.862,       1.33,   1.734,  2.101},
01200     /*19*/{0.861,       1.33,   1.729,  2.093},
01201     /*20*/{0.86,        1.32,   1.725,  2.086},
01202     /*21*/{0.859,       1.32,   1.721,  2.08},
01203     /*22*/{0.858,       1.32,   1.717,  2.074},
01204     /*23*/{0.858,       1.32,   1.714,  2.069},
01205     /*24*/{0.857,       1.32,   1.711,  2.064},
01206     /*25*/{0.856,       1.32,   1.708,  2.06},
01207     /*26*/{0.856,       1,              1.71,   2.06},
01208     /*27*/{0.855,       1,              1.7,    2.05},
01209     /*28*/{0.855,       1,              1.7,    2.05},
01210     /*29*/{0.854,       1,              1.7,    2.04},
01211     /*30*/{0.854,       1,              1.7,    2.04},
01212         
01213               {0.851,   1,              1.68,   2.02},
01214               {0.851,   1,              1.68,   2.02},
01215               {0.851,   1,              1.68,   2.02},
01216               {0.851,   1,              1.68,   2.02},
01217               {0.851,   1,              1.68,   2.02},
01218               {0.851,   1,              1.68,   2.02},
01219               {0.851,   1,              1.68,   2.02},
01220               {0.851,   1,              1.68,   2.02},
01221               {0.851,   1,              1.68,   2.02},
01222         /*40*/{0.851,   1,              1.68,   2.02},
01223                 
01224               {0.848,   1,              1.67,   2},
01225               {0.848,   1,              1.67,   2},
01226               {0.848,   1,              1.67,   2},
01227               {0.848,   1,              1.67,   2},
01228               {0.848,   1,              1.67,   2},
01229               {0.848,   1,              1.67,   2},
01230               {0.848,   1,              1.67,   2},
01231               {0.848,   1,              1.67,   2},
01232               {0.848,   1,              1.67,   2},
01233         /*50*/{0.848,   1,              1.67,   2}
01234         };
01235         
01236     return crtab[alphacat][N-1];
01237 }
01238 
01239 double UCriticalPoint(int N1, int N2, double alpha)
01240 {
01241         static int table1[][10] = {
01242         /*     1   2   3   4   5   6   7   8   9   0*/
01243         /*1*/{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
01244         /*2*/{-1,  4, -1, -1, -1, -1, -1, -1, -1, -1},
01245         /*3*/{-1,  6,  9, -1, -1, -1, -1, -1, -1, -1},
01246         /*4*/{-1,  8, 11, 13, 16, -1, -1, -1, -1, -1},
01247         /*5*/{-1,  9, 13, 16, 20, 23, -1, -1, -1, -1},
01248         /*6*/{-1, 11, 15, 19, 23, 27, -1, -1, -1, -1},
01249         /*7*/{-1, 13, 17, 22, 27, 31, 36, -1, -1, -1},
01250         /*8*/{-1, 14, 19, 25, 30, 35, 40, 45, -1, -1},
01251         /*9*/{ 9, 16, 22, 27, 33, 39, 45, 50, 56, -1},
01252         /*10*/{10,17, 24, 30, 37, 43, 49, 56, 62, 68} 
01253         };
01254     
01255         static int table2[][10] = {
01256         /*     1   2   3   4   5   6   7   8   9   0*/
01257         /*1*/{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
01258         /*2*/{-1, -1, -1, -1, -1, -1, -1, -1, -1, -1},
01259         /*3*/{-1, -1,  9, -1, -1, -1, -1, -1, -1, -1},
01260         /*4*/{-1, -1, 12, 15, 18, -1, -1, -1, -1, -1},
01261         /*5*/{-1, 10, 14, 18, 21, 25, -1, -1, -1, -1},
01262         /*6*/{-1, 12, 16, 21, 25, 29, -1, -1, -1, -1},
01263         /*7*/{-1, 14, 19, 24, 29, 34, 38, -1, -1, -1},
01264         /*8*/{-1, 15, 21, 27, 32, 38, 43, 49, -1, -1},
01265         /*9*/{-1, 17, 23, 30, 36, 42, 48, 54, 60, -1},
01266         /*10*/{-1,19, 26, 33, 39, 46, 53, 60, 66, 73}
01267         };
01268         
01269     int d;
01270     
01271     if (N1 < N2)
01272         {
01273                 if ( alpha <= 0.05 )
01274                 {
01275                 d = table2[N1-1][N2-1];
01276                 }
01277                 else
01278                 {
01279                 d = table1[N1-1][N2-1];
01280                 }
01281         }
01282     else
01283         {
01284                 if ( alpha <= 0.05 )
01285                 {
01286                 d = table2[N2-1][N1-1];
01287                 }
01288                 else
01289                 {
01290                 d = table1[N2-1][N1-1];
01291                 }
01292         }
01293     
01294     if( d == -1 )
01295         {
01296         throw new InvalidArgumentException("UCriticalPoint Invalid sample sizes");
01297         }
01298         return d;
01299 }
01300 
01301 void SampleTest::Rankit(Sample **s1, Sample **s2)
01302 {
01303         List<double> x1;
01304         List<double> x2;
01305         int x;
01306         int N = m_s1.N() + m_s2.N();
01307 
01308         if (N > 50)
01309         {
01310                 throw new Exception("Rankit: Samples too large");
01311         }
01312 
01313         struct _SPair *pairs = new struct _SPair[N];
01314         if ( NULL == pairs )
01315         {
01316                 throw OutOfMemoryException();
01317         }       
01318         for( x = 0; x < N; x++ )
01319         {
01320                 pairs[x].val = -1;
01321                 pairs[x].one = false;
01322         }
01323 
01324         for (x = 0; x < m_s1.N(); x++ )
01325         {
01326                 _ArrayPairInsert (pairs, N, m_s1.Item(x), true);
01327         }
01328 
01329         for (x = 0; x < m_s2.N(); x++ )
01330         {
01331                 _ArrayPairInsert (pairs, N, m_s2.Item(x), true);
01332         }
01333 
01334         // reverse the array
01335         for(x = 0; x < N; x++)
01336         {
01337                 if (x >= N - x)
01338                 {
01339                         break;
01340                 }
01341                 
01342                 double v1 = pairs[x].val;
01343                 bool one1 = pairs[x].one;
01344                 
01345                 pairs[x].val = pairs[(N - 1) - x].val;
01346                 pairs[x].one = pairs[(N - 1) - x].one;
01347                 
01348                 pairs[(N - 1) - x].val = v1;
01349                 pairs[(N - 1) - x].one = one1;
01350         }
01351         ASSERT_MEM(pairs, N * sizeof(struct _SPair));
01352 
01353         // create new packs with the rankit scores as elements
01354         for(x = 0; x < N; x++ )
01355         {
01356                 if (pairs[x].one)
01357                 {
01358                         // my data
01359                         x1.Add( _AvgRankits(pairs, N, x) );
01360                 }
01361                 else
01362                 {
01363                         x2.Add( _AvgRankits(pairs, N, x) );
01364                 }
01365         }
01366         
01367         delete[] pairs;
01368         
01369         if ( NULL == (*s1 = new Sample( &x1 )) )
01370         {
01371                 throw OutOfMemoryException();
01372         }
01373         
01374         if ( NULL == (*s2 = new Sample( &x2 )) )
01375         {
01376                 delete *s1;
01377                 throw OutOfMemoryException();
01378         }
01379 }
01380 
01381 #ifdef DEBUG
01382 void SampleTest::ValidateMem() const
01383 {
01384         m_s1.ValidateMem();
01385         m_s2.ValidateMem();
01386 }
01387 
01388 void SampleTest::CheckMem() const
01389 {
01390         m_s1.CheckMem();
01391         m_s2.CheckMem();
01392 }
01393 #endif