PLearn 0.1
|
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 :