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