PLearn 0.1
|
00001 // -*- C++ -*- 00002 // VecStatsCollector.cc 00003 // 00004 // Copyright (C) 2002 Pascal Vincent 00005 // Copyright (C) 2005 University of Montreal 00006 // Copyright (C) 2007 Xavier Saint-Mleux, ApSTAT Technologies inc. 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: VecStatsCollector.cc 10282 2009-07-21 21:22:48Z plearner $ 00038 ******************************************************* */ 00039 00041 #include "VecStatsCollector.h" 00042 #include "TMat_maths.h" 00043 #include <assert.h> 00044 #include <plearn/base/stringutils.h> 00045 #include <plearn/io/openString.h> 00046 #include <plearn/base/RemoteDeclareMethod.h> 00047 00048 namespace PLearn { 00049 using namespace std; 00050 00051 VecStatsCollector::VecStatsCollector() 00052 : maxnvalues(0), 00053 compute_covariance(false), 00054 epsilon(0.0), 00055 m_window(-1), 00056 m_full_update_frequency(-1), 00057 m_window_nan_code(0), 00058 no_removal_warnings(false), // Window mechanism 00059 sum_non_missing_weights(0), 00060 sum_non_missing_square_weights(0), 00061 m_num_incremental(0) 00062 { } 00063 00064 PLEARN_IMPLEMENT_OBJECT( 00065 VecStatsCollector, 00066 "Collects basic statistics on a vector", 00067 "VecStatsCollector allows to collect statistics on a series of vectors.\n" 00068 "Individual vectors x are presented by calling update(x), and this class will\n" 00069 "collect both individual statistics for each element (as a Vec<StatsCollector>)\n" 00070 "as well as (optionally) compute the covariance matrix." 00071 ); 00072 00073 void VecStatsCollector::declareOptions(OptionList& ol) 00074 { 00075 declareOption( 00076 ol, "maxnvalues", &VecStatsCollector::maxnvalues, 00077 OptionBase::buildoption, 00078 "Maximum number of different values to keep track of for each element.\n" 00079 "If -1, we will keep track of all different values.\n" 00080 "If 0, we will only keep track of global statistics.\n"); 00081 00082 declareOption( 00083 ol, "fieldnames", &VecStatsCollector::fieldnames, OptionBase::buildoption, 00084 "Names of the fields of the vector"); 00085 00086 declareOption( 00087 ol, "compute_covariance", &VecStatsCollector::compute_covariance, OptionBase::buildoption, 00088 "Should we compute and keep the covariance X'X ?"); 00089 00090 declareOption( 00091 ol, "epsilon", &VecStatsCollector::epsilon, OptionBase::buildoption, 00092 "Small regularizing value to be added to the covariance matrix\n" 00093 "estimator, and forwarded to the enclosed vector of StatsCollector.\n" 00094 "This permits dividing by the standard deviation to perform a\n" 00095 "normalization, without fearing a division by zero.\n"); 00096 00097 declareOption( 00098 ol, "window", &VecStatsCollector::m_window, 00099 OptionBase::buildoption, 00100 "If positive, the window restricts the stats computed by this\n" 00101 "VecStatsCollector to the last 'window' observations. This uses the\n" 00102 "VecStatsCollector::remove_observation mechanism.\n" 00103 "Default: -1 (all observations are considered);\n" 00104 " -2 means all observations kept in an ObservationWindow\n"); 00105 00106 declareOption( 00107 ol, "full_update_frequency", &VecStatsCollector::m_full_update_frequency, 00108 OptionBase::buildoption, 00109 "If the window mechanism is used, number of updates at which a full\n" 00110 "update of the underlying StatsCollector is performed. A 'full update'\n" 00111 "is defined as:\n" 00112 "\n" 00113 "- 1. Calling forget()\n" 00114 "- 2. Updating the StatsCollector from all observations in the window.\n" 00115 "\n" 00116 "This is useful for two reasons: 1) when performing a remove-observation\n" 00117 "on a StatsCollector that contains a wide range of values, the\n" 00118 "accumulators for the fourth power may become negative, yielding\n" 00119 "inconsistent estimation. 2) without this option, the statistics\n" 00120 "'FIRST', 'LAST', 'MIN', 'MAX' are not updated properly in the presence\n" 00121 "of a window. To get proper estimation of these statistics, you must\n" 00122 "use the setting 'full_update_frequency=1'.\n" 00123 "\n" 00124 "Default value: -1 (never re-update the StatsCollector from scratch).\n"); 00125 00126 declareOption( 00127 ol, "window_nan_code", &VecStatsCollector::m_window_nan_code, 00128 OptionBase::buildoption, 00129 "How to deal with update vectors containing NaNs with respect to the\n" 00130 "window mechanism.\n" 00131 "\n" 00132 "- 0: Do not check for NaNs (all updates are accounted in the window)\n" 00133 "- 1: If *all* entries of the update vector are NaNs, do not account for\n" 00134 " that observation in the window.\n" 00135 "- 2: If *any* entries of the update vector are NaNs, do not account for\n" 00136 " that observation in the window.\n" 00137 "\n" 00138 " Default: 0" ); 00139 00140 declareOption( 00141 ol, "no_removal_warnings", &VecStatsCollector::no_removal_warnings, 00142 OptionBase::buildoption, 00143 "If the remove_observation mechanism is used (without\n" 00144 "'full_update_frequency=1') and the removed value is equal to one of\n" 00145 "first_, last_, min_ or max_, the default behavior is to warn the user.\n" 00146 "\n" 00147 "To disable this feature, set 'no_removal_warnings' to true.\n" 00148 "\n" 00149 "Default: false (0)." ); 00150 00151 declareOption( 00152 ol, "stats", &VecStatsCollector::stats, OptionBase::learntoption, 00153 "the stats for each element"); 00154 00155 declareOption( 00156 ol, "cov", &VecStatsCollector::cov, OptionBase::learntoption, 00157 "The uncentered and unnormalized covariance matrix (mean not subtracted): X'X"); 00158 00159 declareOption( 00160 ol, "sum_cross", &VecStatsCollector::sum_cross, OptionBase::learntoption, 00161 "Element (i,j) is equal to the (weighted) sum of x_i when both x_i and x_j were observed"); 00162 00163 declareOption( 00164 ol, "sum_cross_weights", &VecStatsCollector::sum_cross_weights, OptionBase::learntoption, 00165 "Element (i,j) is the sum of weights when both x_i and x_j were observed\n" 00166 "(only used when 'compute_covariance' is set to 1)\n"); 00167 00168 declareOption( 00169 ol, "sum_cross_square_weights", &VecStatsCollector::sum_cross_square_weights, OptionBase::learntoption, 00170 "Element (i,j) is the sum of square weights when both x_i and x_j were observed\n" 00171 "(only used when 'compute_covariance' is set to 1)\n"); 00172 00173 declareOption( 00174 ol, "sum_non_missing_weights", &VecStatsCollector::sum_non_missing_weights, OptionBase::learntoption, 00175 "Sum of weights for vectors with no missing value."); 00176 00177 declareOption( 00178 ol, "sum_non_missing_square_weights", &VecStatsCollector::sum_non_missing_square_weights, OptionBase::learntoption, 00179 "Sum of square weights for vectors with no missing value."); 00180 00181 declareOption( 00182 ol, "observation_window", &VecStatsCollector::m_observation_window, 00183 OptionBase::learntoption | OptionBase::nosave | OptionBase::remotetransmit, 00184 "The observation window itself."); 00185 00186 // Now call the parent class' declareOptions 00187 inherited::declareOptions(ol); 00188 } 00189 00191 // declareMethods // 00193 void VecStatsCollector::declareMethods(RemoteMethodMap& rmm) 00194 { 00195 // Insert a backpointer to remote methods; note that this 00196 // different than for declareOptions() 00197 rmm.inherited(inherited::_getRemoteMethodMap_()); 00198 00199 declareMethod( 00200 rmm, "forget", &VecStatsCollector::forget, 00201 (BodyDoc("Clear all previously accumulated statistics.\n"))); 00202 00203 declareMethod( 00204 rmm, "getStat", &VecStatsCollector::getStat, 00205 (BodyDoc("Returns a particular statistic of a particular cost.\n"), 00206 ArgDoc ("statspec", 00207 "A string that is standard statistics specification of the form ex: STAT[fieldname]\n" 00208 "or STAT[fieldnum] where STAT is one of the statistics names understood by\n" 00209 "StatsCollector::getStat. fieldnum start at 0, and fieldnames must have been\n" 00210 "registered with setFieldNames.\n"), 00211 RetDoc ("Requested statistic (a real number)."))); 00212 00213 declareMethod( 00214 rmm, "getMean", &VecStatsCollector::remote_getMean, 00215 (BodyDoc("Return the mean of each field..\n"), 00216 RetDoc ("The vector of means for each field."))); 00217 00218 declareMethod( 00219 rmm, "getVariance", &VecStatsCollector::getVariance, 00220 (BodyDoc("Return the vector of variances of all field..\n"), 00221 RetDoc ("The vector of variance for each field."))); 00222 00223 declareMethod( 00224 rmm, "getStdDev", &VecStatsCollector::getStdDev, 00225 (BodyDoc("Return the vector of standard deviations of all field..\n"), 00226 RetDoc ("The vector of standard deviation for each field."))); 00227 00228 declareMethod( 00229 rmm, "getStdError", &VecStatsCollector::getStdError, 00230 (BodyDoc("Return the vector of standard error of all field..\n"), 00231 RetDoc ("The vector of standard error for each field."))); 00232 00233 declareMethod( 00234 rmm, "getXtX", &VecStatsCollector::getXtX, 00235 (BodyDoc(""), 00236 RetDoc ("Return the matrix XtX "))); 00237 00238 declareMethod( 00239 rmm, "getCovariance", &VecStatsCollector::remote_getCovariance, 00240 (BodyDoc(""), 00241 RetDoc ("Returns the (centered) covariance matrix"))); 00242 00243 declareMethod( 00244 rmm, "getCorrelation", &VecStatsCollector::getCorrelation, 00245 (BodyDoc(""), 00246 RetDoc ("Returns the correlation matrix"))); 00247 00248 declareMethod( 00249 rmm, "setFieldNames", &VecStatsCollector::setFieldNames, 00250 (BodyDoc("Set field names.\n"), 00251 ArgDoc ("fieldnames", 00252 "A vector of strings corresponding to the names of each field" 00253 " in the VecStatsCollector.\n"))); 00254 00255 declareMethod( 00256 rmm, "getFieldNames", &VecStatsCollector::getFieldNames, 00257 (BodyDoc("Get field names.\n"))); 00258 00259 declareMethod( 00260 rmm, "length", &VecStatsCollector::length, 00261 (BodyDoc("Returns the number of statistics collected.\n"), 00262 RetDoc ("=stats.length()"))); 00263 00264 declareMethod( 00265 rmm, "update", &VecStatsCollector::remote_update, 00266 (BodyDoc("Update the stats with gived data.\n"), 00267 ArgDoc ("x"," the new data\n"), 00268 ArgDoc ("weight"," the weight of the data"))); 00269 00270 declareMethod( 00271 rmm, "append", &VecStatsCollector::remote_append, 00272 (BodyDoc("Appends all the StatsCollectors of an " 00273 "existing VecStatsCollector into this one.\n"), 00274 ArgDoc ("vsc","the other VecStatsCollector\n"), 00275 ArgDoc ("fieldname_prefix","prefix concatenated " 00276 "to the existing field names\n"), 00277 ArgDoc ("new_fieldnames","new name for appended fields (overrides prefix)\n"))); 00278 00279 } 00280 00281 int VecStatsCollector::length() const 00282 { 00283 return stats.length(); 00284 } 00285 00286 double VecStatsCollector::getStat(const string& statspec) 00287 { 00288 PStream in = openString(statspec,PStream::plearn_ascii); 00289 string statname; 00290 in.smartReadUntilNext("[", statname); 00291 string fieldname; 00292 in.smartReadUntilNext("]", fieldname); 00293 if(fieldname.empty()) 00294 PLERROR("In VecStatsCollector::getStat - the stat asked is invalid." 00295 "Parsed stat name '%s' with an empty field name.", 00296 statname.c_str()); 00297 int fieldnum = getFieldNum(fieldname); 00298 if(fieldnum<0) 00299 PLERROR("In VecStatsCollector::getStat invalid fieldname: %s;\n" 00300 "valid fieldnames are: %s",fieldname.c_str(), 00301 tostring(fieldnames).c_str()); 00302 00303 // It could be that nothing was accumulated into the stats collector, 00304 // which is different from accessing the "wrong" field. In the first 00305 // case, return MISSING_VALUE 00306 if (stats.length() == 0) 00307 return MISSING_VALUE; 00308 00309 return getStats(fieldnum).getStat(statname); 00310 } 00311 00312 const Mat& VecStatsCollector::getObservations() const 00313 { 00314 PLASSERT( m_window > 0 ); 00315 return m_observation_window->m_observations; 00316 } 00317 00318 const PP<ObservationWindow> 00319 VecStatsCollector::getObservationWindow() const 00320 { 00321 PLASSERT( m_window > 0 ); 00322 return m_observation_window; 00323 } 00324 00325 void VecStatsCollector::setFieldNames(TVec<string> the_fieldnames) 00326 { 00327 fieldnames = the_fieldnames.copy(); 00328 fieldnames_num.clear(); 00329 for (int i=0, n=fieldnames.size() ; i<n ; ++i) 00330 fieldnames_num[fieldnames[i]] = i; 00331 } 00332 00333 int VecStatsCollector::getFieldNum(const string& fieldname_or_num) const 00334 { 00335 map<string,int>::const_iterator it = fieldnames_num.find(fieldname_or_num); 00336 if (it == fieldnames_num.end()) { // not found 00337 if (pl_isnumber(fieldname_or_num)) 00338 return toint(fieldname_or_num); 00339 else 00340 return -1; // unknown field 00341 } 00342 else 00343 return it->second; 00344 } 00345 00346 00348 // update // 00350 void VecStatsCollector::update(const Vec& x, real weight) 00351 { 00352 int n = x.size(); 00353 if(stats.size()==0) 00354 { 00355 stats.resize(n); 00356 for(int k=0; k<n; k++) 00357 { 00358 // TODO It would be cool to have a simple (or automatic) mechanism 00359 // to be able to specify a different value of 'maxnvalues' for each 00360 // StatsCollector (e.g. when only one StatsCollector is meant to 00361 // compute a lift statistics). 00362 stats[k].epsilon = epsilon; 00363 stats[k].maxnvalues = maxnvalues; 00364 stats[k].no_removal_warnings = no_removal_warnings; 00365 stats[k].forget(); 00366 } 00367 if(compute_covariance) 00368 { 00369 cov.resize(n,n); 00370 sum_cross.resize(n,n); 00371 sum_cross_weights.resize(n,n); 00372 sum_cross_square_weights.resize(n,n); 00373 cov.fill(0); 00374 sum_cross.fill(0); 00375 sum_cross_weights.fill(0); 00376 sum_cross_square_weights.fill(0); 00377 } 00378 } 00379 00380 if(stats.size()!=n) 00381 PLERROR("In VecStatsCollector::update - Called update with vector of length " 00382 "%d, while size of stats (and most likely previously seen vector) is " 00383 "%d", n, stats.size()); 00384 00385 // Update the underlying StatsCollectors. If we use the window mechanism 00386 // and we are at a boundary given by m_full_update_frequency, perform a 00387 // full re-update of the StatsCollectors from the saved observations in the 00388 // window. 00389 if ((m_window > 0 || m_window == -2) && m_full_update_frequency > 0 && 00390 ++m_num_incremental >= m_full_update_frequency) 00391 { 00392 for (int k=0 ; k<n ; ++k) 00393 stats[k].forget(); 00394 // Drop oldest observation in window to make room for new observation: 00395 // start at t=1 00396 for (int t=1 ; t<m_observation_window->length() ; ++t) { 00397 Vec obs = m_observation_window->getObs(t); 00398 real w = m_observation_window->getWeight(t); 00399 for (int k=0 ; k<n ; ++k) 00400 stats[k].update(obs[k], w); 00401 } 00402 m_num_incremental = 0; 00403 } 00404 00405 // Incremental update with current observation, as usual 00406 for(int k=0; k<n; ++k) 00407 stats[k].update(x[k], weight); 00408 00409 // Compute covariance if required 00410 if(compute_covariance) { 00411 if (x.hasMissing()) { 00412 // Slower version to handle missing values. 00413 // TODO Could certainly be optimized. 00414 real val_i, val_j; 00415 for (int i = 0; i < n; i++) { 00416 val_i = x[i]; 00417 if (!is_missing(val_i)) { 00418 for (int j = 0; j < n; j++) { 00419 val_j = x[j]; 00420 if (!is_missing(val_j)) { 00421 cov(i,j) += weight * val_i * val_j; 00422 sum_cross(i,j) += weight * val_i; 00423 sum_cross_weights(i,j) += weight; 00424 sum_cross_square_weights(i,j) += weight * weight; 00425 } 00426 } 00427 } 00428 } 00429 } else { 00430 externalProductScaleAcc(cov, x, x, weight); 00431 sum_non_missing_weights += weight; 00432 sum_non_missing_square_weights += weight * weight; 00433 // TODO The two lines below could be optimized with an additional Vec 00434 // storing the sum of weighted x_i for non missing data. 00435 for (int i = 0; i < n; i++) 00436 sum_cross(i) += weight * x[i]; 00437 } 00438 } 00439 00440 // Window mechanism 00441 if ( (m_window > 0 || m_window == -2) && shouldUpdateWindow(x) ) 00442 { 00443 tuple<Vec, real> outdated = m_observation_window->update(x, weight); 00444 Vec& obs = get<0>(outdated); 00445 real w = get<1>(outdated); 00446 00447 // If m_num_incremental==0, we just re-updated the StatsCollectors from 00448 // scratch. In this case, don't call remove_observation. 00449 if(obs.isNotEmpty() && m_window > 0 && 00450 (m_full_update_frequency <= 0 || m_num_incremental > 0)) 00451 { 00452 remove_observation(obs, w); 00453 } 00454 } 00455 } 00456 00457 void VecStatsCollector::remote_update(const Vec& x, real weight) 00458 { 00459 update(x,weight); 00460 } 00461 00462 bool VecStatsCollector::shouldUpdateWindow(const Vec& x) 00463 { 00464 // Avoid dealing with missings if not necessary 00465 if ( m_window_nan_code > 0 ) 00466 { 00467 int count = 0; 00468 Vec::iterator it = x.begin(); 00469 Vec::iterator itend = x.end(); 00470 for(; it!=itend; ++it) 00471 if(is_missing(*it)) 00472 count++; 00473 00474 if ( (m_window_nan_code == 1 && count == x.length()) 00475 || (m_window_nan_code == 2 && count > 0) ) 00476 return false; 00477 } 00478 return true; 00479 } 00480 00482 // remove_observation // 00484 void VecStatsCollector::remove_observation(const Vec& x, real weight) 00485 { 00486 PLASSERT( stats.size() > 0 ); 00487 00488 int n = x.size(); 00489 00490 if(stats.size()!=n) 00491 PLERROR( "In VecStatsCollector: problem, called remove_observation with vector of length %d, " 00492 "while size of stats (and most likeley previously seen vector) is %d", 00493 n, stats.size() ); 00494 00495 for(int k=0; k<n; k++) 00496 { 00497 real obs = x[k]; 00498 stats[k].remove_observation(obs, weight); 00499 //TBA: if ( is_equal(obs, stats[k].min_) ) 00500 //TBA: m_observation_window->columnMin(k, stats[k].min_, stats[k].agemin_); 00501 //TBA: if ( is_equal(obs, stats[k].max_) ) 00502 //TBA: m_observation_window->columnMax(k, stats[k].max_, stats[k].agemax_); 00503 } 00504 00505 // This removes the observation x contribution to the covariance matrix. 00506 if( compute_covariance ) { 00507 if (fast_exact_is_equal(stats[0].nnonmissing(), 0)) { 00508 // We removed the last observation. It may be safer to reset everything 00509 // so that numerical approximations do not lead to negative values for 00510 // statistics that should always be positive. 00511 forget(); 00512 } else { 00513 if (x.hasMissing()) { 00514 // Slower version to handle missing values. 00515 // TODO Could certainly be optimized. 00516 real val_i, val_j; 00517 for (int i = 0; i < n; i++) { 00518 val_i = x[i]; 00519 if (!is_missing(val_i)) { 00520 for (int j = 0; j < n; j++) { 00521 val_j = x[j]; 00522 if (!is_missing(val_j)) { 00523 cov(i,j) -= weight * val_i * val_j; 00524 sum_cross(i,j) -= weight * val_i; 00525 sum_cross_weights(i,j) -= weight; 00526 sum_cross_square_weights(i,j) -= weight * weight; 00527 } 00528 } 00529 } 00530 } 00531 } else { 00532 externalProductScaleAcc(cov, x, x, -weight); 00533 sum_non_missing_weights -= weight; 00534 sum_non_missing_square_weights -= weight * weight; 00535 // TODO The two lines below could be optimized with an additional Vec 00536 // storing the sum of weighted x_i for non missing data. 00537 for (int i = 0; i < n; i++) 00538 sum_cross(i) -= weight * x[i]; 00539 } 00540 } 00541 } 00542 } 00543 00544 00546 void VecStatsCollector::update(const Mat& m) 00547 { 00548 int l = m.length(); 00549 for(int i=0; i<l; i++) 00550 update(m(i)); 00551 } 00552 00554 void VecStatsCollector::update(const Mat& m, const Vec& weights) 00555 { 00556 if (m.length() != weights.size()) 00557 PLERROR("VecStatsCollector::update: matrix height (%d) " 00558 "is incompatible with weights length (%d)", m.length(), 00559 weights.size()); 00560 int l = m.length(); 00561 for(int i=0; i<l; i++) 00562 update(m(i), weights[i]); 00563 } 00564 00565 void VecStatsCollector::build_() 00566 { 00567 if(m_window > 0 || m_window == -2) 00568 { 00569 if ( m_observation_window.isNull() ) 00570 m_observation_window = new ObservationWindow(m_window); 00571 else { 00572 m_observation_window->m_window = m_window; 00573 m_observation_window->forget(); 00574 } 00575 } 00576 00577 if( m_window_nan_code < 0 || m_window_nan_code > 2 ) 00578 PLERROR("The 'window_nan_code' option can only take values 0, 1 or 2."); 00579 } 00580 00581 void VecStatsCollector::build() 00582 { 00583 inherited::build(); 00584 build_(); 00585 } 00586 00587 void VecStatsCollector::forget() 00588 { 00589 stats.resize(0); 00590 cov.resize(0,0); 00591 sum_cross.resize(0,0); 00592 sum_cross_weights.resize(0,0); 00593 sum_cross_square_weights.resize(0,0); 00594 sum_non_missing_weights = 0; 00595 sum_non_missing_square_weights = 0; 00596 00597 // Window mechanism 00598 m_num_incremental = 0; 00599 if (m_window > 0 || m_window == -2) 00600 m_observation_window->forget(); 00601 } 00602 00603 void VecStatsCollector::setWindowSize(int sz) 00604 { 00605 m_window= sz; 00606 if(m_window > 0 || m_window == -2) 00607 { 00608 if(!m_observation_window) 00609 m_observation_window = new ObservationWindow(m_window); 00610 else 00611 { 00612 if(sz == -2) 00613 m_observation_window->unlimited_size= true; 00614 m_observation_window->m_window= sz; 00615 m_observation_window->forget(); 00616 } 00617 } 00618 00619 } 00620 00621 00622 void VecStatsCollector::finalize() 00623 { 00624 int n = stats.size(); 00625 for(int i=0; i<n; i++) 00626 stats[i].finalize(); 00627 } 00628 00630 // getMean // 00632 void VecStatsCollector::getMean(Vec& res) const 00633 { 00634 int n = stats.size(); 00635 res.resize(n); 00636 for(int k=0; k<n; k++) 00637 res[k] = stats[k].mean(); 00638 } 00639 00641 // getVariance // 00643 Vec VecStatsCollector::getVariance() const 00644 { 00645 int n = stats.size(); 00646 Vec res(n); 00647 for(int k=0; k<n; k++) 00648 res[k] = stats[k].variance() + epsilon; 00649 return res; 00650 } 00651 00653 // getStdDev // 00655 Vec VecStatsCollector::getStdDev() const 00656 { 00657 int n = stats.size(); 00658 Vec res(n); 00659 for(int k=0; k<n; k++) 00660 res[k] = sqrt(stats[k].variance() + epsilon); 00661 return res; 00662 } 00663 00665 // getStdError // 00667 Vec VecStatsCollector::getStdError() const 00668 { 00669 int n = stats.size(); 00670 Vec res(n); 00671 for(int k=0; k<n; k++) 00672 res[k] = stats[k].stderror(); 00673 return res; 00674 } 00675 00677 // getCovariance // 00679 void VecStatsCollector::getCovariance(Mat& covar) const { 00680 // Formula used to compute an unbiased estimate of the covariance (note 00681 // that it may not yield a positive semi-definite matrix). 00682 // Notations: 00683 // x(k)_i = i-th coordinate of k-th sample 00684 // sum^i_k f(i,k) = sum over k of f(i,k) for k such that 00685 // x(k)_i is not missing 00686 // sum^i,j_k f(i,j,k) = sum over k of f(i,j,k) for k such that 00687 // neither x(k)_i nor x(k)_j is missing 00688 // w(k) = weight of k-th sample 00689 // cov_i_j = sum^i,j_k w(k) x(k)_i * x(k)_j 00690 // mean_i = (sum^i_k w(k) x(k)_i) / sum^i_k w(k) 00691 // The estimator for element (i,j) of the covariance matrix is then: 00692 // covariance(i,j) = [ cov_i_j + mean_i * mean_j * sum^i,j_k w(k) 00693 // - mean_j * sum^i,j_k w(k) x(k)_i 00694 // - mean_i * sum^i,j_k w(k) x(k)_j 00695 // ] / 00696 // [ sum^i,j_k w(k) 00697 // + (sum^i,j_k w(k) - sum^i_k w(k) - sum^j_k w(k) ) 00698 // * sum^i,j_k w(k)^2 / (sum^i_k w(k) * sum^j_k w(k)) 00699 // ] 00700 // 00701 // If features i and j have never been observed simultaneously, or one of 00702 // the two features has been observed only once, then covariance(i,j) is 00703 // set to MISSING_VALUE. 00704 // The first case occurs when sum^i,j_k w(k) == 0. 00705 // The second case occurs when sum^i,j_k w(k) == sqrt(sum^i,j_k w(k)^2) 00706 // == sum^{i or j}_k w(k) 00707 Vec meanvec; 00708 00709 PLASSERT( compute_covariance && cov.length() == cov.width() ); 00710 int d = cov.length(); 00711 getMean(meanvec); 00712 covar.resize(d,d); 00713 for(int i=0; i<d; i++) { 00714 real sum_weights_i = stats[i].nnonmissing(); 00715 for(int j=i; j<d; j++) { 00716 real sum_weights_j = stats[j].nnonmissing(); 00717 real sum_cross_weights_i_j = sum_cross_weights(i,j) + sum_non_missing_weights; 00718 real sum_cross_square_weights_i_j = sum_cross_square_weights(i,j) 00719 + sum_non_missing_square_weights; 00720 real mean_i = meanvec[i]; 00721 real mean_j = meanvec[j]; 00722 if (fast_exact_is_equal(sum_cross_weights_i_j, 0) || 00723 (fast_is_equal(sum_cross_weights_i_j, 00724 sqrt(sum_cross_square_weights_i_j)) && 00725 (fast_is_equal(sum_cross_weights_i_j, sum_weights_i) || 00726 fast_is_equal(sum_cross_weights_i_j, sum_weights_j)))) 00727 // One of the two cases described above. 00728 covar(i,j) = MISSING_VALUE; 00729 else 00730 covar(i,j) = 00731 (cov(i,j) + mean_i * mean_j * sum_cross_weights_i_j 00732 - mean_j * sum_cross(i,j) 00733 - mean_i * sum_cross(j,i) 00734 ) / 00735 ( sum_cross_weights_i_j 00736 + (sum_cross_weights_i_j - sum_weights_i - sum_weights_j) 00737 * sum_cross_square_weights_i_j 00738 / (sum_weights_i * sum_weights_j) 00739 ); 00740 if (j == i) 00741 covar(i,j) += epsilon; 00742 else 00743 covar(j,i) = covar(i,j); 00744 } 00745 } 00746 } 00747 00748 Mat VecStatsCollector::getCovariance() const 00749 { 00750 Mat covariance; 00751 getCovariance(covariance); 00752 return covariance; 00753 } 00754 00756 // getCorrelation // 00758 Mat VecStatsCollector::getCorrelation() const 00759 { 00760 Mat norm(cov.width(),cov.width()); 00761 externalProduct(norm,getStdDev(),getStdDev()); 00762 return getCovariance()/norm; 00763 } 00764 00765 Vec VecStatsCollector::getAllStats(const string& statname) const 00766 { 00767 const int n = stats.size(); 00768 Vec r(n); 00769 getAllStats(statname, r); 00770 return r; 00771 } 00772 00773 void VecStatsCollector::getAllStats(const string& statname, 00774 Vec& result) const 00775 { 00776 const int n = stats.size(); 00777 result.resize(n); 00778 for (int i=0; i<n; ++i) 00779 result[i] = getStats(i).getStat(statname); 00780 } 00781 00782 void VecStatsCollector::append(const VecStatsCollector& vsc, 00783 const string fieldname_prefix, 00784 const TVec<string>& new_fieldnames) 00785 { 00786 // To avoid problems with fieldnames, ensure we don't start out with too 00787 // many fieldnames, and pad nonexistent fieldnames in *this with "" 00788 fieldnames.resize(stats.size()); 00789 for (int i=fieldnames.size(), n = stats.size() ; i<n ; ++i) 00790 fieldnames.append(""); 00791 00792 stats.append(vsc.stats); 00793 00794 // Take care of field names 00795 if (new_fieldnames.size() > 0) { 00796 PLASSERT( new_fieldnames.size() == vsc.stats.size() ); 00797 fieldnames.append(new_fieldnames); 00798 } 00799 else { 00800 const int n = vsc.stats.size(); 00801 PLASSERT(vsc.fieldnames.size() == n || n == 0); 00802 fieldnames.resize(fieldnames.size(), n); 00803 for (int i=0 ; i<n ; ++i) 00804 fieldnames.append(fieldname_prefix + vsc.fieldnames[i]); 00805 } 00806 setFieldNames(fieldnames); // update map 00807 00808 // Take care of covariance matrix 00809 if (compute_covariance) { 00810 const int oldsize = cov.width(); 00811 const int vscsize = vsc.cov.width(); 00812 PLASSERT( oldsize == cov.length() && vscsize == vsc.cov.length() ); 00813 int new_n = stats.size(); 00814 Mat newcov(new_n, new_n, 0.0); 00815 Mat new_sum_cross(new_n, new_n, 0.0); 00816 Mat new_sum_cross_weights(new_n, new_n, 0.0); 00817 Mat new_sum_cross_square_weights(new_n, new_n, 0.0); 00818 newcov.subMat(0,0,oldsize,oldsize) << cov; 00819 Mat sub = new_sum_cross.subMat(0, 0, oldsize, oldsize); 00820 sub << sum_cross; 00821 sub = new_sum_cross_weights.subMat(0, 0, oldsize, oldsize); 00822 sub << sum_cross_weights; 00823 sub += sum_non_missing_weights; 00824 sum_non_missing_weights = 0; 00825 sub = new_sum_cross_square_weights.subMat(0, 0, oldsize, oldsize); 00826 sub << sum_cross_square_weights; 00827 sub += sum_non_missing_square_weights; 00828 sum_non_missing_square_weights = 0; 00829 if (vsc.compute_covariance) { 00830 newcov.subMat(oldsize,oldsize,vscsize,vscsize) << vsc.cov; 00831 sub = new_sum_cross.subMat(oldsize,oldsize,vscsize,vscsize); 00832 sub << vsc.sum_cross; 00833 sub = new_sum_cross_weights.subMat(oldsize,oldsize,vscsize,vscsize); 00834 sub << vsc.sum_cross_weights; 00835 sub += vsc.sum_non_missing_weights; 00836 sub = new_sum_cross_square_weights.subMat(oldsize,oldsize,vscsize,vscsize); 00837 sub << vsc.sum_cross_square_weights; 00838 sub += vsc.sum_non_missing_square_weights; 00839 } 00840 else { 00841 newcov.subMat(oldsize,oldsize,vscsize,vscsize).fill(MISSING_VALUE); 00842 new_sum_cross.subMat(oldsize, oldsize, vscsize, vscsize).fill(MISSING_VALUE); 00843 new_sum_cross_weights.subMat(oldsize,oldsize,vscsize,vscsize).fill(MISSING_VALUE); 00844 new_sum_cross_square_weights.subMat(oldsize,oldsize,vscsize,vscsize).fill(MISSING_VALUE); 00845 } 00846 cov = newcov; 00847 sum_cross = new_sum_cross; 00848 sum_cross_weights = new_sum_cross_weights; 00849 sum_cross_square_weights = new_sum_cross_square_weights; 00850 } 00851 } 00852 00853 void VecStatsCollector::remote_append(const VecStatsCollector* vsc, 00854 const string fieldname_prefix, 00855 const TVec<string>& new_fieldnames) 00856 { 00857 append(*vsc, fieldname_prefix, new_fieldnames); 00858 } 00859 00860 void VecStatsCollector::merge(VecStatsCollector& other) 00861 { 00862 PLASSERT_MSG(fieldnames == other.fieldnames, 00863 "VecStatsCollector::merge : cannot merge VecStatsCollectors with different fieldnames."); 00864 00865 if(stats.size()==0)//if this one is empty, resize stats before merging 00866 { 00867 int n= other.stats.size(); 00868 stats.resize(n); 00869 for(int k=0; k<n; k++) 00870 { 00871 // TODO It would be cool to have a simple (or automatic) mechanism 00872 // to be able to specify a different value of 'maxnvalues' for each 00873 // StatsCollector (e.g. when only one StatsCollector is meant to 00874 // compute a lift statistics). 00875 stats[k].maxnvalues = maxnvalues; 00876 stats[k].no_removal_warnings = no_removal_warnings; 00877 stats[k].forget(); 00878 } 00879 if(compute_covariance) 00880 { 00881 cov.resize(n,n); 00882 sum_cross.resize(n,n); 00883 sum_cross_weights.resize(n,n); 00884 sum_cross_square_weights.resize(n,n); 00885 cov.fill(0); 00886 sum_cross.fill(0); 00887 sum_cross_weights.fill(0); 00888 sum_cross_square_weights.fill(0); 00889 } 00890 } 00891 00892 PLASSERT_MSG(stats.length() == other.stats.length(), 00893 "VecStatsCollector::merge : cannot merge VecStatsCollectors with different stats length."); 00894 00895 if(m_window != 0 || other.m_window != 0) 00896 { 00897 PLASSERT_MSG(m_window != 0 && other.m_window != 0, 00898 "VecStatsCollector::merge : either none or both should have an observation window"); 00899 PP<ObservationWindow> oow= other.m_observation_window; 00900 for(int i= 0; i < oow->length(); ++i) 00901 update(oow->getObs(i), oow->getWeight(i)); 00902 return; // avoid extra indentation 00903 } 00904 00905 for(int i= 0; i < stats.length(); ++i) 00906 stats[i].merge(other.stats[i]); 00907 00908 if(compute_covariance) 00909 { 00910 for(int i= 0; i < cov.length(); ++i) 00911 for(int j= 0; j < cov.width(); ++j) 00912 { 00913 cov(i,j)+= other.cov(i,j); 00914 sum_cross(i,j)+= other.sum_cross(i,j); 00915 sum_cross_weights(i,j)+= other.sum_cross_weights(i,j); 00916 sum_cross_square_weights(i,j)+= other.sum_cross_square_weights(i,j); 00917 } 00918 sum_non_missing_weights+= other.sum_non_missing_weights; 00919 sum_non_missing_square_weights+= other.sum_non_missing_square_weights; 00920 } 00921 00922 } 00923 00924 void VecStatsCollector::makeDeepCopyFromShallowCopy(CopiesMap& copies) 00925 { 00926 inherited::makeDeepCopyFromShallowCopy(copies); 00927 deepCopyField(fieldnames, copies); 00928 deepCopyField(stats, copies); 00929 deepCopyField(cov, copies); 00930 deepCopyField(sum_cross, copies); 00931 deepCopyField(sum_cross_weights, copies); 00932 deepCopyField(sum_cross_square_weights, copies); 00933 deepCopyField(m_observation_window, copies); 00934 } 00935 00936 } // end of namespace PLearn 00937 00938 00939 /* 00940 Local Variables: 00941 mode:c++ 00942 c-basic-offset:4 00943 c-file-style:"stroustrup" 00944 c-file-offsets:((innamespace . 0)(inline-open . 0)) 00945 indent-tabs-mode:nil 00946 fill-column:79 00947 End: 00948 */ 00949 // vim: filetype=cpp:expandtab:shiftwidth=4:tabstop=8:softtabstop=4:encoding=utf-8:textwidth=79 :