PLearn 0.1
VMat_basic_stats.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // VMat_basic_stats.cc
00004 //
00005 // Copyright (C) 2004 Pascal Vincent
00006 // Copyright (C) 2005 University of Montreal
00007 //
00008 // Redistribution and use in source and binary forms, with or without
00009 // modification, are permitted provided that the following conditions are met:
00010 //
00011 //  1. Redistributions of source code must retain the above copyright
00012 //     notice, this list of conditions and the following disclaimer.
00013 //
00014 //  2. Redistributions in binary form must reproduce the above copyright
00015 //     notice, this list of conditions and the following disclaimer in the
00016 //     documentation and/or other materials provided with the distribution.
00017 //
00018 //  3. The name of the authors may not be used to endorse or promote
00019 //     products derived from this software without specific prior written
00020 //     permission.
00021 //
00022 // THIS SOFTWARE IS PROVIDED BY THE AUTHORS ``AS IS'' AND ANY EXPRESS OR
00023 // IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES
00024 // OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN
00025 // NO EVENT SHALL THE AUTHORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
00026 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED
00027 // TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00028 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00029 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00030 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00031 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00032 //
00033 // This file is part of the PLearn library. For more information on the PLearn
00034 // library, go to the PLearn Web site at www.plearn.org
00035 
00036 /* *******************************************************
00037  * $Id: VMat_basic_stats.cc 8916 2008-04-30 14:12:24Z nouiz $
00038  ******************************************************* */
00039 
00040 // Authors: Pascal Vincent
00041 
00044 #include <plearn/base/Object.h>
00045 #include "VMat_basic_stats.h"
00046 //#include "VMat.h"
00047 #include "MemoryVMatrix.h"
00048 #include "ShiftAndRescaleVMatrix.h"
00049 //#include <plearn/math/TMat_maths.h>
00050 #include <plearn/math/stats_utils.h>
00051 #include <plearn/math/VecStatsCollector.h>
00052 #include <plearn/math/TMat_maths.h>
00053 //#include <plearn/sys/PLMPI.h>
00054 
00055 namespace PLearn {
00056 using namespace std;
00057 
00059 // computeWeightedMean //
00061 void computeWeightedMean(const Vec& weights, const VMat& d, Vec& meanvec)
00062 {
00063     VecStatsCollector sc;
00064     int n = d->length();
00065     if (weights.length() != n)
00066         PLERROR("In computeWeightedMean - weights.length() != d->length()");
00067     Vec row(d->width());
00068     for (int i = 0; i < n; i++) {
00069         d->getRow(i, row);
00070         sc.update(row, weights[i]);
00071     }
00072     sc.getMean(meanvec);
00073 }
00074 
00076 // computeRange //
00078 void computeRange(const VMat& d, Vec& minvec, Vec& maxvec)
00079 {
00080     int n = d->length();
00081     int w = d->width();
00082     minvec.resize(w);
00083     maxvec.resize(w);
00084     VecStatsCollector sc;
00085     Vec row(w);
00086     for (int i = 0; i < n; i++) {
00087         d->getRow(i, row);
00088         sc.update(row);
00089     }
00090     for (int j = 0; j < w; j++) {
00091         minvec[j] = sc.getStats(j).min();
00092         maxvec[j] = sc.getStats(j).max();
00093     }
00094 }
00095 
00097 // computeRowMean //
00099 void computeRowMean(const VMat& d, Vec& meanvec)
00100 {
00101     int n = d->length();
00102     meanvec.resize(n);
00103     Vec samplevec(d->width());
00104     for(int i = 0; i < n; i++)
00105     {
00106         d->getRow(i,samplevec);
00107         meanvec[i] = mean(samplevec);
00108     }
00109 }
00110 
00112 // computeMean //
00114 void computeMean(const VMat& d, Vec& meanvec)
00115 {
00116     Vec constant_weight(d->length(), 1.0);
00117     computeWeightedMean(constant_weight, d, meanvec);
00118 }
00119 
00121 // computeBasicStats //
00123 Mat computeBasicStats(const VMat& m)
00124 {
00125     // TODO Use StatsCollector instead ?
00126     Vec v(m.width());
00127     real* vdata = v.data();
00128     Mat stats(10,m.width());
00129     Vec mean_row = stats(MEAN_ROW);
00130     Vec stddev_row = stats(STDDEV_ROW);
00131     Vec min_row = stats(MIN_ROW);
00132     Vec max_row = stats(MAX_ROW);
00133     Vec nmissing_row = stats(NMISSING_ROW);
00134     Vec nzero_row = stats(NZERO_ROW);
00135     Vec npositive_row = stats(NPOSITIVE_ROW);
00136     Vec nnegative_row = stats(NNEGATIVE_ROW);
00137     Vec meanpos_row = stats(MEANPOS_ROW);
00138     Vec stddevpos_row = stats(STDDEVPOS_ROW);
00139     min_row.fill(FLT_MAX);
00140     max_row.fill(-FLT_MAX);
00141 
00142     for(int i=0; i<m.length(); i++)
00143     {
00144         m->getRow(i,v);
00145         for(int j=0; j<v.length(); j++)
00146         {
00147             real val = vdata[j];
00148             if(is_missing(val))
00149                 nmissing_row[j]++;
00150             else
00151             {
00152                 if(val<min_row[j])
00153                     min_row[j] = val;
00154                 if(val>max_row[j])
00155                     max_row[j] = val;
00156 
00157                 if(fast_exact_is_equal(val, 0.))
00158                     nzero_row[j]++;
00159                 else if(val>0.)
00160                 {
00161                     npositive_row[j]++;
00162                     mean_row[j] += val;
00163                     stddev_row[j] += val*val;
00164                     meanpos_row[j] += val;
00165                     stddevpos_row[j] += val*val;
00166                 }
00167                 else // val < 0.
00168                 {
00169                     nnegative_row[j]++;
00170                     mean_row[j] += val;
00171                     stddev_row[j] += val*val;
00172                 }
00173             }
00174         }
00175     }
00176     for(int j=0; j<stats.width(); j++)
00177     {
00178         real nnonmissing = nzero_row[j]+nnegative_row[j]+npositive_row[j];
00179         mean_row[j] /= nnonmissing;
00180         meanpos_row[j] /= npositive_row[j];
00181         stddev_row[j] = sqrt(stddev_row[j]/nnonmissing - square(mean_row[j]));
00182         stddevpos_row[j] = sqrt(stddevpos_row[j]/npositive_row[j] - square(meanpos_row[j]));
00183     }
00184     return stats;
00185 }
00186 
00188 // computeConditionalMeans //
00190 // mean = sum/n
00191 // variance = (sumsquare-square(sum)/n)/(n-1)
00192 // stddev_of_mean = sqrt(variance/n);
00193 // mse = sumsquare/n - square(sum/n)
00194 // stddev_of_mse = variance*sqrt(2./n);
00195 TVec<Mat> computeConditionalMeans(const VMat& trainset, int targetsize, Mat& basic_stats)
00196 {
00197     if(!basic_stats)
00198         basic_stats = computeBasicStats(trainset);
00199 
00200     int inputsize = trainset.width()-targetsize;
00201     TVec<Mat> a(inputsize);
00202     for(int j=0; j<inputsize; j++)
00203     {
00204         real minval = basic_stats(MIN_ROW,j);
00205         real maxval = basic_stats(MAX_ROW,j);
00206         if(is_integer(minval) && is_integer(maxval) && maxval-minval<400)
00207         {
00208             a[j] = Mat(int(maxval-minval+1),2+targetsize*4);
00209             for(int k=0; k<a[j].length(); k++)
00210                 a[j](k,0) = minval+k;
00211         }
00212     }
00213 
00214     Vec row(trainset.width());
00215     Vec input = row.subVec(0,inputsize);
00216     Vec target = row.subVec(inputsize,targetsize);
00217     for(int i=0; i<trainset.length(); i++)
00218     {
00219         trainset->getRow(i,row);
00220         for(int j=0; j<inputsize; j++)
00221         {
00222             Mat& m = a[j];
00223             if(m.isNotEmpty())
00224             {
00225                 int k = int(input[j]-basic_stats(MIN_ROW,j));
00226                 Vec m_k = m(k);
00227                 m_k[1]++;
00228                 for(int l=0; l<targetsize; l++)
00229                 {
00230                     real targetval = target[l];
00231                     m_k[2+4*l] += targetval;
00232                     m_k[3+4*l] += square(targetval);
00233                 }
00234             }
00235         }
00236     }
00237 
00238     // postprocessing:
00239     for(int j=0; j<inputsize; j++)
00240     {
00241         Mat& m = a[j];
00242         if(m.isNotEmpty())
00243         {
00244             for(int k=0; k<m.length(); k++)
00245             {
00246                 Vec m_k = m(k);
00247                 real n = m_k[1];
00248                 if(n>0.)
00249                 {
00250                     // replace sum by mean and sumsquare by variance
00251                     for(int l=0; l<targetsize; l++)
00252                     {
00253                         real sum = m_k[2+4*l];
00254                         real sumsquare = m_k[3+4*l];
00255                         real mean = sum/n;
00256                         real variance = (sumsquare-square(sum)/n)/(n-1);
00257                         real mse = sumsquare/n - square(sum/n);
00258                         real stddev_of_mean = sqrt(variance/n);
00259                         real stddev_of_mse = variance*sqrt(2./n);
00260                         m_k[2+4*l] = mean;
00261                         m_k[3+4*l] = stddev_of_mean;
00262                         m_k[4+4*l] = mse;
00263                         m_k[5+4*l] = stddev_of_mse;
00264                     }
00265                 }
00266             }
00267         }
00268     }
00269 
00270     return a;
00271 }
00272 
00274 // computeMeanAndVariance //
00276 void computeMeanAndVariance(const VMat& d, Vec& meanvec, Vec& variancevec,
00277                             real epsilon)
00278 {
00279     VecStatsCollector sc;
00280     sc.epsilon = epsilon;
00281     sc.build();
00282     int n = d->length();
00283     Vec row(d->width());
00284     for (int i = 0; i < n; i++) {
00285         d->getRow(i, row);
00286         sc.update(row);
00287     }
00288     sc.getMean(meanvec);
00289     variancevec.resize(d->width());
00290     variancevec << sc.getVariance();
00291 }
00292 
00294 // computeInputMean //
00296 void computeInputMean(const VMat& d, Vec& meanvec)
00297 {
00298     VecStatsCollector sc;
00299     int n = d->length();
00300     Vec input, target;
00301     real weight;
00302     for (int i = 0; i < n; i++) {
00303         d->getExample(i, input, target, weight);
00304         sc.update(input, weight);
00305     }
00306     sc.getMean(meanvec);
00307 }
00308 
00310 // computeInputMeanAndCovar //
00312 void computeInputMeanAndCovar(const VMat& d, Vec& meanvec, Mat& covarmat,
00313                               real epsilon)
00314 {
00315     PLASSERT( d->inputsize() >= 0 );
00316     VecStatsCollector sc;
00317     sc.compute_covariance = true;
00318     sc.epsilon = epsilon;
00319     sc.build();
00320     int n = d->length();
00321     Vec input, target;
00322     real weight;
00323     for (int i = 0; i < n; i++) {
00324         d->getExample(i, input, target, weight);
00325         sc.update(input, weight);
00326     }
00327     sc.getMean(meanvec);
00328     sc.getCovariance(covarmat);
00329 }
00330 
00332 // computeInputMeanAndVariance //
00334 void computeInputMeanAndVariance(const VMat& d, Vec& meanvec, Vec& var,
00335                                  real epsilon)
00336 {
00337     PLASSERT( d->inputsize() >= 0 );
00338     VecStatsCollector sc;
00339     sc.epsilon=epsilon;
00340     sc.build();
00341     int n = d->length();
00342     Vec input, target;
00343     real weight;
00344     for (int i = 0; i < n; i++) {
00345         d->getExample(i, input, target, weight);
00346         sc.update(input, weight);
00347     }
00348     sc.getMean(meanvec);
00349     var.resize(d->inputsize());
00350     var << sc.getVariance();
00351 }
00352 
00354 // computeInputMeanAndStddev //
00356 void computeInputMeanAndStddev(const VMat& d, Vec& meanvec, Vec& stddev,
00357                                real epsilon)
00358 {
00359     computeInputMeanAndVariance(d, meanvec, stddev, epsilon);
00360     for (int i = 0; i < stddev.length(); i++) {
00361 #ifdef BOUNDCHECK
00362         if (stddev[i] < 0)
00363             PLERROR("In computeInputMeanAndStddev - The computed variance should be >= 0");
00364 #endif
00365         stddev[i] = sqrt(stddev[i]);
00366     }
00367 }
00368 
00370 // computeWeightedMeanAndCovar //
00372 void computeWeightedMeanAndCovar(const Vec& weights, const VMat& d, Vec& meanvec, Mat& covarmat,
00373                                  real epsilon)
00374 {
00375     VecStatsCollector sc;
00376     sc.compute_covariance = true;
00377     sc.epsilon = epsilon;
00378     sc.build();
00379     int n = d->length();
00380     Vec row(d->width());
00381     for (int i = 0; i < n; i++) {
00382         d->getRow(i, row);
00383         sc.update(row, weights[i]);
00384     }
00385     sc.getMean(meanvec);
00386     sc.getCovariance(covarmat);
00387 }
00388 
00390 // computeMeanAndCovar //
00392 void computeMeanAndCovar(const VMat& d, Vec& meanvec, Mat& covarmat, real epsilon)
00393 {
00394     VecStatsCollector sc;
00395     sc.compute_covariance = true;
00396     sc.epsilon = epsilon;
00397     sc.build();
00398     int n = d->length();
00399     Vec row(d->width());
00400     for (int i = 0; i < n; i++) {
00401         d->getRow(i, row);
00402         sc.update(row);
00403     }
00404     sc.getMean(meanvec);
00405     sc.getCovariance(covarmat);
00406 
00407     /* Commented out old code that had an optimized MPI version, but was probably
00408        not used anymore.
00409 
00410        int w = m->width();
00411        int l = m->length();
00412        meanvec.resize(w);
00413        covarmat.resize(w,w);
00414 
00415        MemoryVMatrix* memvm = dynamic_cast<MemoryVMatrix*>((VMatrix*)m);
00416        if(memvm)
00417        computeMeanAndCovar(m->toMat(), meanvec, covarmat);
00418        else
00419        {
00420        meanvec.clear();
00421        covarmat.clear();
00422        Vec v(w);
00423 
00424        ProgressBar progbar("Computing covariance",l);
00425 
00426        if(USING_MPI && PLMPI::synchronized && PLMPI::size>1)
00427        { //!<  Parallel implementation
00428        #if USING_MPI
00429        PLMPI::synchronized = false;
00430 
00431        if(!covarmat.isCompact())
00432        PLERROR("In computeMeanAndCovar: MPI implementation cannot handle non-compact covariance matrices, please pass a compact matrix");
00433 
00434        // temporary storages for mpi
00435        Vec meanvec_b(meanvec.length());
00436        Mat covarmat_b(covarmat.length(),covarmat.width());
00437 
00438        for(int i=PLMPI::rank; i<l; i+=PLMPI::size)
00439        {
00440        m->getRow(i,v);
00441        meanvec_b += v;
00442        externalProductAcc(covarmat_b, v, v);
00443        progbar(i);
00444        }
00445 
00446        MPI_Reduce(meanvec_b.data(), meanvec.data(), meanvec.length(), PLMPI_REAL, MPI_SUM, 0, MPI_COMM_WORLD);
00447        MPI_Bcast(meanvec.data(), meanvec.length(), PLMPI_REAL, 0, MPI_COMM_WORLD);
00448        MPI_Reduce(covarmat_b.data(), covarmat.data(), covarmat.size(), PLMPI_REAL, MPI_SUM, 0, MPI_COMM_WORLD);
00449        MPI_Bcast(covarmat.data(), covarmat.size(), PLMPI_REAL, 0, MPI_COMM_WORLD);
00450 
00451        PLMPI::synchronized = true;
00452        #endif
00453        }
00454        else //!<  default sequential implementation
00455        {
00456        for(int i=0; i<l; i++)
00457        {
00458        m->getRow(i,v);
00459        meanvec += v;
00460        externalProductAcc(covarmat, v, v);
00461        progbar(i);
00462        }
00463        }
00464 
00465        // get the real averages and covariances, and priors
00466        meanvec /= real(l);
00467        covarmat /= real(l);
00468        externalProductScaleAcc(covarmat,meanvec,meanvec,real(-1.));
00469        }
00470     */
00471 }
00472 
00473 void computeCovar(const VMat& d, const Vec& mu, Mat& covarmat, real epsilon)
00474 {
00475     int w = d->width();
00476     int l = d->length();
00477     covarmat.resize(w,w);
00478     covarmat.clear();
00479     Vec samplevec(w);
00480     Vec diffvec(w);
00481     Mat sqdiffmat(w,w);
00482     for(int i=0; i<l; i++)
00483     {
00484         d->getRow(i,samplevec);
00485         samplevec -= mu;
00486         externalProductAcc(covarmat, samplevec, samplevec);
00487     }
00488     covarmat /= l-1;
00489     addToDiagonal(covarmat, epsilon);
00490 }
00491 
00492 void computeInputCovar(const VMat& d, const Vec& mu, Mat& covarmat, real epsilon)
00493 {
00494     PLASSERT( d->inputsize() >= 0 );
00495     int w = d->inputsize();
00496     int l = d->length();
00497     covarmat.resize(w,w);
00498     covarmat.clear();
00499     Vec input(w);
00500     Vec target;
00501     real weight;
00502     Vec diffvec(w);
00503     Mat sqdiffmat(w,w);
00504     real weightsum = 0;
00505     for(int i=0; i<l; i++)
00506     {
00507         d->getExample(i, input, target, weight);
00508         input -= mu;
00509         externalProductScaleAcc(covarmat, input, input, weight);
00510         weightsum += weight;
00511     }
00512     covarmat *= real(1./weightsum);
00513     addToDiagonal(covarmat, epsilon);
00514 }
00515 
00516 
00518 // computeMeanAndStddev //
00520 void computeMeanAndStddev(const VMat& d, Vec& meanvec, Vec& stddevvec,
00521                           real epsilon)
00522 {
00523     computeMeanAndVariance(d, meanvec, stddevvec, epsilon);
00524     for(int i=0; i<stddevvec.length(); i++)
00525         stddevvec[i] = sqrt(stddevvec[i]);
00526 }
00527 
00529 // autocorrelation_function //
00531 void autocorrelation_function(const VMat& data, Mat& acf)
00532 {
00533     int T = data.length();
00534     int N = data.width();
00535     acf.resize(T-2, N);
00536 
00537     for(int delta=0; delta < T-2; delta++)
00538     {
00539         Vec sumT(N);
00540         Vec sumD(N);
00541         TVec<Vec> products(N);
00542 
00543         // t = delta
00544         for(int k=0; k < N; k++)
00545         {
00546             real ts = data(delta, k);
00547             real ds = data(0, k);
00548 
00549             sumT[k] = ts;
00550             sumD[k] = ds;
00551 
00552             products[k].resize(3);
00553             products[k][0] = ts*ts;
00554             products[k][1] = ds*ds;
00555             products[k][2] = ts*ds;
00556         }
00557 
00558         for(int t=delta+1; t < T; t++)
00559         {
00560             for(int k=0; k < N; k++)
00561             {
00562                 real ts = data(t, k);
00563                 real ds = data(t-delta, k);
00564 
00565                 sumT[k] += ts;
00566                 sumD[k] += ds;
00567 
00568                 products[k][0] += ts*ts;
00569                 products[k][1] += ds*ds;
00570                 products[k][2] += ts*ds;
00571             }
00572         }
00573 
00574         // Actual computation of the correlation
00575         for(int k=0; k < N; k++)
00576         {
00577             int count = T-delta;
00578             real multiplied_var_t = products[k][0] - square(sumT[k])/count;
00579             real multiplied_var_d = products[k][1] - square(sumD[k])/count;
00580             acf(delta, k) = (products[k][2] - sumT[k]*sumD[k]/count) / sqrt(multiplied_var_t * multiplied_var_d);
00581         }
00582     }
00583 }
00584 
00585 
00587 // normalize //
00589 VMat normalize(const VMat& d, const Vec& meanvec, const Vec& stddevvec)
00590 {
00591     int inputsize = meanvec.length();
00592 
00593     Vec shiftvec(d.width(), 0.0);
00594     shiftvec.subVec(0,inputsize) << meanvec;
00595     negateElements(shiftvec);
00596 
00597     Vec scalevec(d.width(), 1.0);
00598     scalevec.subVec(0,inputsize) << stddevvec;
00599     invertElements(scalevec);
00600 
00601     return new ShiftAndRescaleVMatrix(d, shiftvec, scalevec);
00602 }
00603 
00605 // normalize //
00607 VMat normalize(const VMat& d, int inputsize, int ntrain)
00608 {
00609     Vec meanvec(inputsize);
00610     Vec stddevvec(inputsize);
00611     computeMeanAndStddev(d.subMat(0,0,ntrain,inputsize), meanvec, stddevvec);
00612     return normalize(d, meanvec, stddevvec);
00613 }
00614 
00616 // normalize //
00618 VMat normalize(VMat d, int inputsize)
00619 {
00620     return normalize(d, inputsize, d.length());
00621 }
00622 
00624 // correlations //
00626 void correlations(const VMat& x, const VMat& y, Mat& r, Mat& pvalues, bool ignore_missing)
00627 {
00628     TMat<int> n_nonmissing; // Store the number of non-missing values for each pair.
00629     int n=x.length();
00630     if (n!=y.length())
00631         PLERROR("correlations: x and y must have the same length");
00632     int wx=x.width();
00633     int wy=y.width();
00634     r.resize(wx,wy);
00635     r.clear();
00636     Mat sxy(wx,wy);
00637     Vec sx2(wx);
00638     Vec sy2(wy);
00639     Vec sx(wx);
00640     Vec sy(wy);
00641     Vec xt(wx);
00642     Vec yt(wy);
00643     Mat sy_m, sx_m, sy2_m, sx2_m;
00644     if (ignore_missing) {
00645         n_nonmissing.resize(wx, wy);
00646         sy_m.resize(wx, wy);
00647         sy2_m.resize(wx, wy);
00648         sx_m.resize(wx, wy);
00649         sx2_m.resize(wx, wy);
00650         n_nonmissing.fill(0);
00651         sy_m.fill(0);
00652         sy2_m.fill(0);
00653         sx_m.fill(0);
00654         sx2_m.fill(0);
00655     }
00656     for (int t=0;t<n;t++)
00657     {
00658         x->getRow(t,xt);
00659         y->getRow(t,yt);
00660         for (int j=0;j<wy;j++)
00661         {
00662             real ytj = yt[j];
00663             if (!ignore_missing) {
00664 #ifdef BOUNDCHECK
00665                 if (is_missing(ytj))
00666                     PLWARNING("In correlations - You should not compute correlations "
00667                               "with missing values and 'ignore_ missing' set to false");
00668 #endif
00669                 sy[j] += ytj;
00670                 sy2[j] += ytj*ytj;
00671             }
00672             for (int i=0;i<wx;i++)
00673             {
00674                 real xti = xt[i];
00675                 if (ignore_missing) {
00676                     if (!is_missing(ytj) && !is_missing(xti)) {
00677                         sy_m(i,j) += ytj;
00678                         sy2_m(i,j) += ytj * ytj;
00679                         sx_m(i,j) += xti;
00680                         sx2_m(i,j) += xti * xti;
00681                         sxy(i,j) += xti * ytj;
00682                         n_nonmissing(i,j)++;
00683                     }
00684                 } else {
00685 #ifdef BOUNDCHECK
00686                     if (is_missing(xti))
00687                         PLWARNING("In correlations - You should not compute correlations "
00688                                   "with missing values and 'ignore_ missing' set to false");
00689 #endif
00690                     sxy(i,j) += xti*ytj;
00691                     sx[i] += xti;
00692                     sx2[i] += xti*xti;
00693                 }
00694             }
00695         }
00696     }
00697     for (int i=0;i<wx;i++)
00698         for (int j=0;j<wy;j++)
00699         {
00700             real nv; // = n * variance of x
00701             if (ignore_missing) {
00702                 nv = sx2_m(i,j) - sx_m(i,j) / real(n_nonmissing(i,j)) * sx_m(i,j);
00703             } else {
00704                 nv = sx2[i] - sx[i]/real(n)*sx[i];
00705             }
00706             if (nv>0) // don't bother if variance is 0
00707                 if (ignore_missing)
00708                     r(i,j) = (n_nonmissing(i,j)*sxy(i,j)-sx_m(i,j)*sy_m(i,j)) /
00709                         sqrt( (n_nonmissing(i,j)*sx2_m(i,j)-sx_m(i,j)*sx_m(i,j)) *
00710                               (n_nonmissing(i,j)*sy2_m(i,j)-sy_m(i,j)*sy_m(i,j)));
00711                 else
00712                     r(i,j) = (n*sxy(i,j)-sx[i]*sy[j])/sqrt((n*sx2[i]-sx[i]*sx[i])*(n*sy2[j]-sy[j]*sy[j]));
00713             else
00714                 r(i,j) = 0;
00715             if (r(i,j)<-1.01 || r(i,j)>1.01)
00716                 PLWARNING("correlation: weird correlation coefficient, %f for %d-th input, %d-target",
00717                           r(i,j),i,j);
00718         }
00719     pvalues.resize(wx, wy);
00720     for (int i=0;i<wx;i++)
00721         for (int j=0;j<wy;j++)
00722             if (ignore_missing)
00723                 pvalues(i,j) = testNoCorrelationAsymptotically(r(i,j),n_nonmissing(i,j));
00724             else
00725                 pvalues(i,j) = testNoCorrelationAsymptotically(r(i,j),n);
00726 }
00727 
00728 } // end of namespace PLearn
00729 
00730 
00731 /*
00732   Local Variables:
00733   mode:c++
00734   c-basic-offset:4
00735   c-file-style:"stroustrup"
00736   c-file-offsets:((innamespace . 0)(inline-open . 0))
00737   indent-tabs-mode:nil
00738   fill-column:79
00739   End:
00740 */
00741 // 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