PLearn 0.1
stats_utils.cc
Go to the documentation of this file.
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 :
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines