PLearn 0.1
|
00001 // -*- C++ -*- 00002 00003 // PLearn (A C++ Machine Learning Library) 00004 // Copyright (C) 2003 Pascal Vincent, Yoshua Bengio 00005 00006 // Redistribution and use in source and binary forms, with or without 00007 // modification, are permitted provided that the following conditions are met: 00008 // 00009 // 1. Redistributions of source code must retain the above copyright 00010 // notice, this list of conditions and the following disclaimer. 00011 // 00012 // 2. Redistributions in binary form must reproduce the above copyright 00013 // notice, this list of conditions and the following disclaimer in the 00014 // documentation and/or other materials provided with the distribution. 00015 // 00016 // 3. The name of the authors may not be used to endorse or promote 00017 // products derived from this software without specific prior written 00018 // permission. 00019 // 00020 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR 00021 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES 00022 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN 00023 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 00024 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED 00025 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR 00026 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF 00027 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING 00028 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS 00029 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 00030 // 00031 // This file is part of the PLearn library. For more information on the PLearn 00032 // library, go to the PLearn Web site at www.plearn.org 00033 00034 00035 /* ******************************************************* 00036 * $Id: stats_utils.cc 9789 2008-12-16 21:35:14Z nouiz $ 00037 * This file is part of the PLearn library. 00038 ******************************************************* */ 00039 00040 00043 #include "stats_utils.h" 00044 #include "TMat_maths.h" 00045 #include "pl_erf.h" 00046 #include "random.h" 00047 #include <plearn/base/RemoteDeclareMethod.h> 00048 00049 namespace PLearn { 00050 using namespace std; 00051 00053 // SpearmanRankCorrelation // 00055 TMat<int> SpearmanRankCorrelation(const VMat &x, const VMat& y, Mat& r, bool ignore_missing) 00056 { 00057 TMat<int> result; 00058 int n=x.length(); 00059 if (n!=y.length()) 00060 PLERROR("SpearmanRankCorrelation: x and y must have the same length"); 00061 int wx=x.width(); 00062 int wy=y.width(); 00063 r.resize(wx,wy); 00064 r.clear(); 00065 Mat y_ranks; 00066 TVec<int> n_ynonmissing = computeRanks(y.toMat(),y_ranks, ignore_missing); 00067 Mat x_rank(n,1); 00068 //real rank_normalization = sqrt(1.0/(n*n-1.0)); 00069 real rank_normalization = 12.0/(n*(n-1.0)*(n-2.0)); 00070 real half = n*0.5; 00071 // Vectors used only when 'ignore_missing' is true. 00072 Vec rank_normalization_miss; 00073 Vec half_miss; 00074 Mat y_copy, y_rankj, xi_placeholder; 00075 if (ignore_missing) { 00076 rank_normalization_miss.resize(wy); 00077 half_miss.resize(wy); 00078 result.resize(wx, wy); 00079 y_copy.resize(n, 1); 00080 y_rankj.resize(n, 1); 00081 xi_placeholder.resize(x.length(), 1); 00082 } 00083 ProgressBar pb("Computing Spearman rank correlation", wx); 00084 for (int i=0;i<wx; pb.update(++i)) 00085 { 00086 Mat xi = x.column(i).toMat(); 00087 const Vec& r_i = r(i); 00088 00089 if (ignore_missing) { 00090 // Ignoring missing values is more complex, because we need to ignore them 00091 // both in x and y. 00092 for (int j = 0; j < wy; j++) { 00093 // We replace values in x_i associated to missing values in y by missing 00094 // values, so that they are not taken into account when computing ranks. 00095 // We need to make a copy to ensure we do not destroy any data. 00096 if (n_ynonmissing[j] < y.length()) { 00097 Mat xi_back = xi; // Backup of original xi. 00098 xi = xi_placeholder; 00099 xi << xi_back; 00100 for (int k = 0; k < y_ranks.length(); k++) 00101 if (fast_exact_is_equal(y_ranks(k,j), -1)) // -1 rank <-> missing value 00102 xi(k,0) = MISSING_VALUE; 00103 } 00104 TVec<int> n_nonmissing = computeRanks(xi, x_rank, ignore_missing); 00105 n = n_nonmissing[0]; 00106 result(i,j) = n; 00107 rank_normalization = 12.0/(n*(n-1.0)*(n-2.0)); 00108 half = n*0.5; 00109 y_rankj << y_ranks.column(j); 00110 // We need to recompute y's ranks if there were missing values in x. 00111 if (n < n_ynonmissing[j]) { 00112 y_copy << y.column(j).toMat(); 00113 for (int k = 0; k < x_rank.length(); k++) 00114 if (fast_exact_is_equal(x_rank(k,0), -1)) 00115 y_copy(k,0) = MISSING_VALUE; 00116 computeRanks(y_copy, y_rankj, ignore_missing); 00117 } else 00118 y_rankj << y_ranks.column(j); 00119 for (int k = 0; k < x_rank.length(); k++) { 00120 real x_r = x_rank(k,0); 00121 real y_r = y_rankj(k,0); 00122 if (!fast_exact_is_equal(x_r, -1)) // -1 rank <-> missing value 00123 r_i[j] += (x_r - half) * (y_r - half) * rank_normalization; 00124 } 00125 } 00126 } else { 00127 // Compute the rank of the i-th column of x 00128 computeRanks(xi,x_rank, ignore_missing); 00129 // Compute the Spearman rank correlation coefficient: 00130 for (int k=0;k<n;k++) 00131 for (int j=0;j<wy;j++) 00132 { 00133 //real delta = (x_rank(k,0) - y_ranks(k,j))*rank_normalization; 00134 // r_i[j] += delta*delta; 00135 r_i[j] += (x_rank(k,0) - half) * (y_ranks(k,j)-half) * rank_normalization; 00136 } 00137 } 00138 for (int j=0;j<wy;j++) 00139 if (r_i[j]<-1.01 || r_i[j]>1.01) 00140 PLWARNING("SpearmanRankCorrelation: weird correlation coefficient, %f for %d-th input, %d-target", 00141 r_i[j],i,j); 00142 } 00143 return result; 00144 } 00145 00151 real testNoCorrelationAsymptotically(real r, int n) 00152 { 00153 real fz = fabs(r)*sqrt(n-1.0); 00154 return (1-gauss_01_cum(fz)) + gauss_01_cum(-fz); 00155 } 00156 00157 void testSpearmanRankCorrelationPValues(const VMat &x, const VMat& y, Mat& pvalues, bool ignore_missing) 00158 { 00159 Mat r; 00160 testSpearmanRankCorrelation(x,y,r,pvalues, ignore_missing); 00161 } 00162 00163 void testSpearmanRankCorrelation(const VMat &x, const VMat& y, Mat& r, Mat& pvalues, bool ignore_missing) 00164 { 00165 int n = x.length(); 00166 TMat<int> n_mat = SpearmanRankCorrelation(x,y,r, ignore_missing); 00167 pvalues.resize(r.length(),r.width()); 00168 for (int i=0;i<r.length();i++) 00169 for (int j=0;j<r.width();j++) 00170 pvalues(i,j) = testNoCorrelationAsymptotically(r(i,j), ignore_missing ? n_mat(i,j) : n); 00171 } 00172 00173 00175 real max_cdf_diff(Vec& v1, Vec& v2) 00176 { 00177 int n1 = v1.length(); 00178 int n2 = v2.length(); 00179 real inv_n1 = 1./n1; 00180 real inv_n2 = 1./n2; 00181 sortElements(v1); 00182 sortElements(v2); 00183 int i1=0; 00184 int i2=0; 00185 real maxdiff = 0; 00186 00187 if(n1==0 && n2==0) 00188 { 00189 PLWARNING("In max_cdf_diff(Vec, Vec) - both vector are empty!"); 00190 return 0; 00191 } 00192 else if (n1==0 || n2==0) 00193 return 1; 00194 00195 for(;;) 00196 { 00197 00198 if(v1[i1]<v2[i2]) 00199 { 00200 i1++; 00201 if(i1+1==n1) 00202 break; 00203 } 00204 else if(fast_exact_is_equal(v1[i1],v2[i2])) 00205 { 00206 i1++;i2++; 00207 if(i2==n2) 00208 break; 00209 else if(i1+1==n1) 00210 break; 00211 continue; 00212 } 00213 else 00214 { 00215 i2++; 00216 if(i2==n2) 00217 break; 00218 } 00219 00220 if ((i1>0 && fast_exact_is_equal(v1[i1], v1[i1-1])) || 00221 (i2>0 && fast_exact_is_equal(v2[i2], v2[i2-1])) || 00222 (v1[i1]<v2[i2] && v1[i1+1]<v2[i2])) 00223 continue; // to deal with discrete-valued variables: only look at "changing-value" places 00224 00225 real F1 = inv_n1*i1; 00226 real F2 = inv_n2*i2; 00227 real diff = fabs(F1-F2); 00228 if(diff>maxdiff) 00229 maxdiff = diff; 00230 00231 // perr << "v1[" << i1 << "]=" << v1[i1] << "; v2[" << i2 << "]=" << v2[i2] << "; F1=" << F1 << "; F2=" << F2 << "; diff=" << diff << endl; 00232 } 00233 00234 return maxdiff; 00235 } 00236 00237 00238 /*************************************************************/ 00239 /* 00240 Return the probability that the Kolmogorov-Smirnov statistic 00241 D takes the observed value or greater, given the null hypothesis 00242 that the distributions that are compared are 00243 really identical. N is the effective number of samples 00244 used for comparing the distributions. 00245 The argument conv gives the precision with which 00246 this probability is computed. A value above 10 does not bring 00247 much improvement. Note that the statistic D can 00248 be obtained as follows: 00249 00250 Comparing two empirical distributions from data sets D_1 and D_2: 00251 Let F_1(x) the empirical cumulative distribution of D_1 of size N_1, and 00252 let F_2(x) the empirical cumulative distribution of D_2 of size N_2. Then 00253 00254 D = max_x | F_1(x) - F_2(x) | 00255 00256 and the effective N is N_1 N_2 / (N_1 + N_2). 00257 00258 Comparing a theoretical distribution F and a data set D of size N with 00259 empirical cumulative distribution F_N: 00260 00261 D = max_x | F(x) - F_N(x) | 00262 00263 This function returns the following 00264 00265 P(D > observed d | same distributions) estimated by 00266 2 sum_{k=1}^{infty} (-1)^{k-1} exp(-2k^2 a^2) 00267 00268 where a = sqrt(D*(sqrt(N)+0.12+0.11/sqrt(N))) 00269 00270 Ref: Stephens, M.A. (1970), Journal of the Royal Statistical Society B, vol. 32, pp. 115-122. 00271 00272 */ 00273 real KS_test(real D, real N, int conv) 00274 { 00275 int k; 00276 real res = 0.0; 00277 real sn = sqrt((double)N); 00278 real ks = D*(sn+0.12+0.11/sn); 00279 real ks2 = ks*ks; 00280 for (k=1;k<=conv;k++) { 00281 real x = ((k % 2) ? 1 : -1) * exp( -2 * ks2 * k * k ); 00282 if (k==conv) 00283 res += 0.5*x; 00284 else 00285 res += x; 00286 } 00287 return 2 * res; 00288 } 00289 00290 void KS_test(Vec& v1, Vec& v2, int conv, real& D, real& p_value) 00291 { 00292 int n1 = v1.length(); 00293 int n2 = v2.length(); 00294 real N = (n1/real(n1+n2))*n2; 00295 D = max_cdf_diff(v1, v2); 00296 p_value = KS_test(D,N,conv); 00297 } 00298 00300 void KS_test(const VMat& m1, const VMat& m2, const int conv, Vec& Ds, Vec& p_values, 00301 const bool report_progress) 00302 { 00303 m1->compatibleSizeError(m2); 00304 Ds.resize(m1->width()); 00305 p_values.resize(m1->width()); 00306 PP<ProgressBar> pbar; 00307 if (report_progress) 00308 pbar = new ProgressBar("Computing Kologorov Smirnow two sample test", 00309 m1->width()); 00310 00311 #pragma omp parallel default(none) shared(pbar) 00312 { 00313 Vec row1(m1->length()); 00314 Vec row2(m2->length()); 00315 #pragma omp for 00316 for(int col = 0;col<m1->width();col++) 00317 { 00318 row1->resize(m1->length()); 00319 row2->resize(m2->length()); 00320 #pragma omp critical 00321 m1->getColumn(col,row1);//not threadsafe! 00322 #pragma omp critical 00323 m2->getColumn(col,row2);//not threadsafe! 00324 remove_missing_inplace(row1); 00325 remove_missing_inplace(row2); 00326 real D; 00327 real p_value; 00328 KS_test(row1,row2,conv,D,p_value); 00329 Ds[col]=D; 00330 p_values[col]=p_value; 00331 if (report_progress) 00332 pbar->update(col); 00333 } 00334 } 00335 } 00336 real KS_test(Vec& v1, Vec& v2, int conv) 00337 { 00338 real D, ks_stat; 00339 KS_test(v1,v2,conv,D, ks_stat); 00340 return ks_stat; 00341 } 00342 00343 tuple<real,real> remote_KS_test(Vec& v1, Vec& v2, int conv) 00344 { 00345 real D, pvalue; 00346 KS_test(v1,v2,conv,D, pvalue); 00347 return make_tuple(D, pvalue); 00348 } 00349 00350 tuple<Vec,Vec> remote_KS_tests(VMat& m1, VMat& m2, int conv) 00351 { 00352 Vec Ds, pvalues; 00353 KS_test(m1, m2, conv, Ds, pvalues); 00354 return make_tuple(Ds, pvalues); 00355 } 00356 00357 real paired_t_test(Vec u, Vec v) 00358 { 00359 int n = u.length(); 00360 if( v.length() != n ) 00361 { 00362 PLWARNING("paired_t_test: " 00363 "Can't make a paired t-test on to unequally lengthed vectors (%d != %d).", 00364 n, v.length()); 00365 return MISSING_VALUE; 00366 } 00367 00368 real ubar = mean(u); 00369 real vbar = mean(v); 00370 Vec u2 = u - ubar; 00371 Vec v2 = v - vbar; 00372 00373 return (ubar - vbar) * sqrt( n*(n-1) / sumsquare(u2-v2)); 00374 } 00375 00384 void DirichletEstimatorMMoments(const Mat& p, Vec& alpha) 00385 { 00386 static Vec mean_p, mean_p2, var_p; // NON-REENTRANT CODE 00387 int N=p.width(); 00388 alpha.resize(N); 00389 mean_p.resize(N); 00390 mean_p2.resize(N); 00391 var_p.resize(N); 00392 columnMean(p, mean_p); 00393 columnSumOfSquares(p, mean_p2); 00394 mean_p2 *= real(1.0/N); 00395 columnVariance(p, var_p, mean_p); 00396 real log_sum_alpha = 0; 00397 for (int i=0;i<N;i++) 00398 log_sum_alpha += safeflog(mean_p[i]*(1-mean_p[i])/var_p[i]-1); 00399 log_sum_alpha /= (N-1); 00400 multiply(mean_p, exp(log_sum_alpha), alpha); 00401 } 00402 00415 void DirichletEstimatorMaxLik(const Mat& p, Vec alpha) 00416 { 00417 DirichletEstimatorMMoments(p,alpha); 00418 // int N=alpha.length(); Commented out to remove compiler warning. 00419 PLERROR("In DirichletEstimatorMaxLik - Not implemented yet"); 00420 // Have a look at Tom Minka's paper on estimating Dirichlet parameters... 00421 // TO BE IMPLEMENTED 00422 } 00423 00424 00425 00426 00427 BEGIN_DECLARE_REMOTE_FUNCTIONS 00428 00429 declareFunction("KS_test", &remote_KS_test, 00430 (BodyDoc("Returns result of Kolmogorov-Smirnov test between 2 samples.\n"), 00431 ArgDoc ("v1","Vec1: first distr."), 00432 ArgDoc ("v2","Vec2: second distr."), 00433 ArgDoc ("conv","precision"), 00434 RetDoc ("tuple of (D, p-value)"))); 00435 00436 declareFunction("KS_tests", &remote_KS_tests, 00437 (BodyDoc("Returns result of Kolmogorov-Smirnov test between 2 VMats, for each column.\n"), 00438 ArgDoc ("m1","VMat1: first distr."), 00439 ArgDoc ("m2","VMat2: second distr."), 00440 ArgDoc ("conv","precision"), 00441 RetDoc ("tuple of (Ds, p-values)"))); 00442 00443 END_DECLARE_REMOTE_FUNCTIONS 00444 00445 00446 00447 } // end of namespace PLearn 00448 00449 /* 00450 00451 // Test code... 00452 00453 #include "random.h" 00454 #include <plearn/display/Gnuplot.h> 00455 00456 using namespace PLearn; 00457 00458 // should plot a uniform distribution 00459 void verify_ks(int n1=1000, int n2=1000, int k=100) 00460 { 00461 Vec v1(n1); 00462 Vec v2(n2); 00463 Vec ks(k); 00464 00465 for(int i=0; i<k; i++) 00466 { 00467 fill_random_normal(v1, 0, 1); 00468 fill_random_normal(v2, 0.1, 1); 00469 // fill_random_uniform(v2, -0.5, 0.5); 00470 00471 ks[i] = KS_test(v1,v2); 00472 perr << '.'; 00473 } 00474 00475 Gnuplot gp; 00476 gp.plotcdf(ks); 00477 char s[100]; 00478 cin.getline(s, 100); 00479 } 00480 00481 int main() 00482 { 00483 Vec v1(5); 00484 v1 << "1 2 5 9 14"; 00485 00486 Vec v2(6); 00487 v2 << "-1 4 12 14 25 3"; 00488 00489 real md = max_cdf_diff(v1,v2); 00490 00491 int n1 = v1.length(); 00492 int n2 = v2.length(); 00493 00494 pout << md << endl; 00495 pout << KS_test(md, n1*n2/real(n1+n2)) << endl; 00496 00497 verify_ks(); 00498 00499 return 0; 00500 } 00501 00502 */ 00503 00504 00505 00506 /* 00507 Local Variables: 00508 mode:c++ 00509 c-basic-offset:4 00510 c-file-style:"stroustrup" 00511 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00512 indent-tabs-mode:nil 00513 fill-column:79 00514 End: 00515 */ 00516 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :