PLearn 0.1
PDistribution.cc
Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // PDistribution.cc
00004 //
00005 // Copyright (C) 2003 Pascal Vincent
00006 // Copyright (C) 2004-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: PDistribution.cc 9789 2008-12-16 21:35:14Z nouiz $
00038  ******************************************************* */
00039 
00042 #include "PDistribution.h"
00043 #include <plearn/base/tostring.h>
00044 #include <plearn/math/TMat_maths.h>
00045 #include <plearn/math/PRandom.h>
00046 
00047 namespace PLearn {
00048 using namespace std;
00049 
00051 // PDistribution //
00053 PDistribution::PDistribution():
00054     delta_curve(0.1),
00055     predictor_size(0),
00056     predicted_size(-1),
00057     n_predictor(-1),
00058     n_predicted(-1),
00059     lower_bound(0.),
00060     upper_bound(0.),
00061     n_curve_points(-1),
00062     outputs_def("l")
00063 {
00064     random_gen = new PRandom();
00065 }
00066 
00067 PLEARN_IMPLEMENT_OBJECT(
00068     PDistribution,
00069     "Base class for PLearn probability distributions.",
00070     "PDistributions derive from PLearner, as some of them may be fitted to data by\n"
00071     "training, but they have additional methods allowing e.g. to compute density\n"
00072     "or generate data points.\n"
00073     "\n"
00074     "By default, a PDistribution may be conditional to a predictor part x, in\n"
00075     "order to represent the conditional distribution of P(Y | X = x). An\n"
00076     "unconditional distribution should derive from UnconditionalDistribution as it\n"
00077     "has a simpler interface.\n"
00078     "\n"
00079     "Since we want to be able to compute for instance P(Y = y | X = x), both the\n"
00080     "predictor part 'x' and the predicted part 'y' must be considered as input\n"
00081     "from the PLearner framework point of view. Thus one must specify the size of\n"
00082     "the predictor part by the 'predictor_size' option, and the size of the\n"
00083     "predicted by the 'predicted_size' option, satisfying the following equality:\n"
00084     "\n"
00085     "    predictor_size + predicted_size == inputsize  (1)\n"
00086     "\n"
00087     "Optionally, 'predictor_size' or 'predicted_size' (but not both) may be set to\n"
00088     "-1, and the PDistribution will automatically guess the other size so that\n"
00089     "equation (1) is satisfied (actually, in order to preserve the user-provided\n"
00090     "values of 'predictor_size' and 'predicted_size', the guessed values are\n"
00091     "stored in the learnt options 'n_predictor' and 'n_predicted'). This way,\n"
00092     "unconditional distributions can be created by setting 'predictor_size' to 0\n"
00093     "and 'predicted_size' to -1.\n"
00094     "\n"
00095     "The default implementations of the learner-type methods for computing\n"
00096     "outputs and costs work as follows:\n"
00097     "  - the 'outputs_def' option allows to choose what is in the output\n"
00098     "    (e.g. log density, expectation, ...)\n"
00099     "  - the cost is a vector of size 1 containing only the negative log-\n"
00100     "    likelihood (NLL), i.e. -log(P(y|x)).\n"
00101     "\n"
00102     "For conditional distributions, the input must always be made of both the\n"
00103     "'predictor' part (x) and the 'predicted' part (y), even if the output may not\n"
00104     "need the predicted part (e.g. to compute E[Y | X = x]).  The exception is\n"
00105     "when computeOutput(..) needs to be called successively with the same value of\n"
00106     "'x': in this case, after a first call with both 'x' and 'y', one may only\n"
00107     "provide 'y' as input in later calls, and 'x' will be assumed to be\n"
00108     "unchanged. Or, alternatively, one can set the 'predictor_part' option first,\n"
00109     "either through the options system or using the setPredictor(..) method.\n"
00110     );
00111 
00113 // declareOptions //
00115 void PDistribution::declareOptions(OptionList& ol)
00116 {
00117 
00118     // Build options.
00119 
00120     declareOption(
00121         ol, "outputs_def", &PDistribution::outputs_def,
00122                            OptionBase::buildoption,
00123         "Defines what will be given in output. This is a string where the\n"
00124         "characters have the following meaning:\n"
00125         "- 'l' : log_density\n"
00126         "- 'd' : density\n"
00127         "- 'c' : cdf\n"
00128         "- 's' : survival_fn\n"
00129         "- 'e' : expectation\n"
00130         "- 'v' : variance.\n"
00131         "\n"
00132         "If these options are specified in lower case they give the value\n"
00133         "associated with a given observation. In upper case, a curve is\n"
00134         "evaluated at regular intervals and produced in output (as a\n"
00135         "histogram). For 'L', 'D', 'C', 'S', it is the predicted part that\n"
00136         "varies, while for 'E' and 'V' it is the predictor part (for\n"
00137         "conditional distributions).\n"
00138         "The number of curve points is given by the 'n_curve_points' option.\n"
00139         "Note that the upper case letters only work for scalar variables, in\n"
00140         "order to produce a one-dimensional curve."
00141         );
00142     // TODO Make it TVec<string> for better clarity?
00143 
00144     declareOption(ol, "predictor_size",  &PDistribution::predictor_size,
00145                                   OptionBase::buildoption,
00146         "The (user-provided) size of the predictor x in p(y|x). A value of\n"
00147         "-1 means the algorithm should find it out by itself.");
00148 
00149     declareOption(ol, "predicted_size", &PDistribution::predicted_size,
00150                                         OptionBase::buildoption,
00151         "The (user-provided) size of the predicted y in p(y|x). A value of\n"
00152         "-1 means the algorithm should find it out by itself.");
00153 
00154     declareOption(ol, "predictor_part", &PDistribution::predictor_part,
00155                                         OptionBase::buildoption,
00156         "In conditional distributions, the predictor part (x in P(Y|X=x)).\n");
00157 
00158     declareOption(ol, "n_curve_points", &PDistribution::n_curve_points,
00159                                         OptionBase::buildoption,
00160         "The number of points for which the output is evaluated when\n"
00161         "outputs_defs is upper case (produces a histogram).\n"
00162         "The lower_bound and upper_bound options specify where the curve\n"
00163         "begins and ends. Note that these options (upper case letters) only\n"
00164         "work for scalar variables.");
00165 
00166     declareOption(ol, "lower_bound",  &PDistribution::lower_bound,
00167                                       OptionBase::buildoption,
00168         "The lower bound of scalar Y values to compute a histogram of the\n"
00169         "distribution when upper case outputs_def are specified.");
00170 
00171     declareOption(ol, "upper_bound",  &PDistribution::upper_bound,
00172                                       OptionBase::buildoption,
00173         "The upper bound of scalar Y values to compute a histogram of the\n"
00174         "distribution when upper case outputs_def are specified.");
00175 
00176     // Learnt options.
00177 
00178     declareOption(ol, "n_predictor",  &PDistribution::n_predictor,
00179                                       OptionBase::learntoption,
00180         "The (true) size of the predictor x in p(y|x). If 'predictor_size'\n"
00181         "is non-negative, 'n_predictor' is set to 'predictor_size'.\n"
00182         "Otherwise, it is set to the data dimension minus 'predicted_size'.");
00183 
00184     declareOption(ol, "n_predicted",  &PDistribution::n_predicted,
00185                                       OptionBase::learntoption,
00186         "The (true) size of the predicted y in p(y|x). If 'predicted_size'\n"
00187         "is non-negative, 'n_predicted' is set to 'predicted_size'.\n"
00188         "Otherwise, it is set to the data dimension minus 'predictor_size'.");
00189 
00190     // Now call the parent class' declareOptions
00191     inherited::declareOptions(ol);
00192 
00193 }
00194 
00196 // declareMethods //
00198 void PDistribution::declareMethods(RemoteMethodMap& rmm)
00199 {
00200     // Insert a backpointer to remote methods; note that this
00201     // different from declareOptions()
00202     rmm.inherited(inherited::_getRemoteMethodMap_());
00203 
00204     declareMethod(
00205         rmm, "log_density", &PDistribution::log_density,
00206         (BodyDoc("Compute the log density of a data point"),
00207          ArgDoc("sample", "The data point"),
00208          RetDoc("The log density")));
00209 
00210     declareMethod(
00211         rmm, "generate", &PDistribution::remote_generate,
00212         (BodyDoc("Generate a sample"),
00213          RetDoc("The generated sample")));
00214 }
00215 
00217 // build //
00219 void PDistribution::build()
00220 {
00221     inherited::build();
00222     build_();
00223 }
00224 
00226 // build_ //
00228 void PDistribution::build_()
00229 {
00230     // Reset the random number generator seed.
00231     resetGenerator(seed_);
00232 
00233     // Typical code for a PDistribution: the class makes the operations it
00234     // needs when the predictor and predicted sizes are defined, and when the
00235     // predictor is defined. In the build_() method, it should not call the
00236     // parent's methods since they should have already been called during the
00237     // parent's build.
00238     PDistribution::setPredictorPredictedSizes(predictor_size, predicted_size,
00239                                               false);
00240     PDistribution::setPredictor(predictor_part, false);
00241 
00242     // Set the step between two points in the output curve.
00243     if (n_curve_points > 0)
00244         delta_curve = (upper_bound - lower_bound) / real(n_curve_points);
00245 }
00246 
00248 // computeOutput //
00250 void PDistribution::computeOutput(const Vec& input, Vec& output) const
00251 {
00252     // TODO Add an output to generate samples.
00253 
00254     // Set the 'predictor' (x in P(Y = y| X=x)) and 'predicted' (y) parts.
00255     splitCond(input);
00256 
00257     string::size_type l = outputs_def.length();
00258     output.resize(outputsize());
00259 
00260     int k = 0;
00261     for(unsigned int i=0; i<l; i++)
00262     {
00263         switch(outputs_def[i])
00264         {
00265         case 'l':
00266             output[k++] = log_density(predicted_part);
00267             break;
00268         case 'd':
00269             output[k++] = density(predicted_part);
00270             break;
00271         case 'c':
00272             output[k++] = cdf(predicted_part);
00273             break;
00274         case 's':
00275             output[k++] = survival_fn(predicted_part);
00276             break;
00277         case 'e':
00278             store_expect = output.subVec(k, n_predicted);
00279             expectation(store_expect);
00280             k += n_predicted;
00281             break;
00282         case 'v':
00283             store_cov =
00284                 output.subVec(k, square(n_predicted)).toMat(n_predicted,n_predicted);
00285             variance(store_cov);
00286             k += square(n_predicted);
00287             break;
00288         case 'E':
00289         case 'V':
00290             if (n_predicted > 1)
00291                 PLERROR("In PDistribution::computeOutput - Can only plot "
00292                         "histogram of expectation or variance for "
00293                         "one-dimensional expected part");
00294             if (n_predicted == 0)
00295                 PLERROR("In PDistribution::computeOutput - Cannot plot "
00296                         "histogram of expectation or variance for "
00297                         "unconditional distributions");
00298         case 'L':
00299         case 'D':
00300         case 'C':
00301         case 'S':
00302             real t;
00303             store_result.resize(1);
00304             store_result[0] = lower_bound;
00305             for (int j = 0; j < n_curve_points; j++) {
00306                 switch(outputs_def[i]) {
00307                 case 'L':
00308                     t = log_density(store_result);
00309                     break;
00310                 case 'D':
00311                     t = density(store_result);
00312                     break;
00313                 case 'C':
00314                     t = cdf(store_result);
00315                     break;
00316                 case 'S':
00317                     t = survival_fn(store_result);
00318                     break;
00319                 case 'E':
00320                     setPredictor(store_result);
00321                     expectation(store_expect);
00322                     t = store_expect[0];
00323                     break;
00324                 case 'V':
00325                     setPredictor(store_result);
00326                     store_cov = store_expect.toMat(1,1);
00327                     variance(store_cov);
00328                     t = store_expect[0];
00329                     break;
00330                 default:
00331                     PLERROR("In PDistribution::computeOutput - This should "
00332                             "never happen");
00333                     t = 0; // To make the compiler happy.
00334                 }
00335                 output[j + k] = t;
00336                 store_result[0] += delta_curve;
00337             }
00338             k += n_curve_points;
00339             break;
00340         default:
00341             // Maybe a subclass knows about this output?
00342             // TODO This is quite ugly. See how to do this better.
00343             unknownOutput(outputs_def[i], input, output, k);
00344             break;
00345         }
00346     }
00347 }
00348 
00350 // computeCostsFromOutputs //
00352 void PDistribution::computeCostsFromOutputs(const Vec& input, const Vec& output,
00353                                             const Vec& target, Vec& costs) const
00354 {
00355     costs.resize(1);
00356     char c = outputs_def[0];
00357     if(c == 'l')
00358     {
00359         costs[0] = -output[0];
00360     }
00361     else if(c == 'd')
00362     {
00363         costs[0] = -pl_log(output[0]);
00364     }
00365     else
00366         PLERROR("In PDistribution::computeCostsFromOutputs currently can only "
00367                 "compute' NLL cost from log likelihood or density returned as "
00368                 "first output");
00369 }
00370 
00372 // getTestCostNames //
00374 TVec<string> PDistribution::getTestCostNames() const
00375 {
00376     TVec<string> nll_cost;
00377     if (nll_cost.isEmpty())
00378         nll_cost.append("NLL");
00379     return nll_cost;
00380 }
00381 
00383 // getTrainCostNames //
00385 TVec<string> PDistribution::getTrainCostNames() const
00386 {
00387     // Default = no train cost computed. This may be overridden in subclasses.
00388     TVec<string> no_cost;
00389     return no_cost;
00390 }
00391 
00393 // generateN //
00395 void PDistribution::generateN(const Mat& Y) const
00396 {
00397     Vec v;
00398     if (Y.width() != n_predicted)
00399         PLERROR("In PDistribution::generateN - Matrix width (%d) differs from "
00400                 "n_predicted (%d)", Y.width(), n_predicted);
00401     int N = Y.length();
00402     PP<ProgressBar> pb =
00403         report_progress ? new ProgressBar("Generating samples", N)
00404                         : NULL;
00405     for(int i=0; i<N; i++)
00406     {
00407         v = Y(i);
00408         generate(v);
00409         if (pb)
00410             pb->update(i);
00411     }
00412 }
00413 
00415 // makeDeepCopyFromShallowCopy //
00417 void PDistribution::makeDeepCopyFromShallowCopy(CopiesMap& copies)
00418 {
00419     inherited::makeDeepCopyFromShallowCopy(copies);
00420     deepCopyField(store_expect,       copies);
00421     deepCopyField(store_result,       copies);
00422     deepCopyField(store_cov,          copies);
00423     deepCopyField(predictor_part,     copies);
00424     deepCopyField(predicted_part,     copies);
00425 }
00426 
00428 // outputsize //
00430 int PDistribution::outputsize() const
00431 {
00432     int l = 0;
00433     for (size_t i=0; i<outputs_def.length(); i++) {
00434         if (outputs_def[i]=='L' || outputs_def[i]=='D' || outputs_def[i]=='C'
00435          || outputs_def[i]=='S' || outputs_def[i]=='E' || outputs_def[i]=='V')
00436             l+=n_curve_points;
00437         else if (outputs_def[i]=='e')
00438             l += n_predicted;
00439         else if (outputs_def[i]=='v')
00440             // Variance is full (n x n) matrix.
00441             l += square(n_predicted);
00442         else l++;
00443     }
00444     return l;
00445 }
00446 
00447 
00449 // remote_generate //
00451 Vec PDistribution::remote_generate()
00452 {
00453     Vec sample;
00454     generate(sample);
00455     return sample;
00456 }
00457 
00459 // resetGenerator //
00461 void PDistribution::resetGenerator(long g_seed)
00462 {
00463     if (g_seed != 0) {
00464         seed_ = g_seed;
00465         random_gen->manual_seed(g_seed);
00466     }
00467 }
00468 
00470 // setPredictor //
00472 void PDistribution::setPredictor(const Vec& predictor, bool call_parent) const
00473 {
00474     // Default behavior: only fill 'predictor_part' with first elements of
00475     // 'predictor'.
00476     PLASSERT( predictor.length()      >= n_predictor );
00477     PLASSERT( predictor_part.length() == n_predictor );
00478     if (predictor != predictor_part)
00479         predictor_part << predictor.subVec(0, n_predictor);
00480 }
00481 
00483 // setPredictorPredictedSizes //
00485 bool PDistribution::setPredictorPredictedSizes(int the_predictor_size,
00486                                                int the_predicted_size,
00487                                                bool call_parent)
00488 {
00489     PLASSERT( (the_predictor_size  >= 0 || the_predictor_size  == -1) &&
00490             (the_predicted_size >= 0 || the_predicted_size == -1) );
00491     int backup_n_predictor = n_predictor;
00492     int backup_n_predicted = n_predicted;
00493     n_predictor = predictor_size = the_predictor_size;
00494     n_predicted = predicted_size = the_predicted_size;
00495     if (n_predictor < 0) {
00496         if (n_predicted < 0)
00497             PLERROR("In PDistribution::setPredictorPredictedSizes - You need"
00498                     "to specify at least one non-negative value");
00499         if (inputsize_ >= 0) {
00500             if (n_predicted > inputsize_)
00501                 PLERROR("In PDistribution::setPredictorPredictedSizes - "
00502                         "'n_predicted' (%d) cannot be > inputsize (%d)",
00503                         n_predicted, inputsize_);
00504             n_predictor = inputsize_ - n_predicted;
00505         }
00506     } else if (n_predicted < 0) {
00507         if (inputsize_ >= 0) {
00508             if (n_predictor > inputsize_)
00509                 PLERROR("In PDistribution::setPredictorPredictedSizes - "
00510                         "'n_predictor' (%d) cannot be > inputsize (%d)",
00511                         n_predictor, inputsize_);
00512             n_predicted = inputsize_ - n_predictor;
00513         }
00514     }
00515     if (inputsize_ >= 0 && n_predictor + n_predicted != inputsize_)
00516         PLERROR("In PDistribution::setPredictorPredictedSizes - n_predictor "
00517                 "(%d) + n_predicted (%d) != inputsize (%d)",
00518                 n_predictor, n_predicted, inputsize_);
00519     if (n_predictor >= 0)
00520         predictor_part.resize(n_predictor);
00521     if (n_predicted >= 0)
00522         predicted_part.resize(n_predicted);
00523     if (!call_parent)
00524         return false;
00525     else
00526         return (n_predictor != backup_n_predictor ||
00527                 n_predicted != backup_n_predicted);
00528 }
00529 
00531 // splitCond //
00533 void PDistribution::splitCond(const Vec& input) const {
00534     if (n_predictor == 0 || (n_predictor > 0 && input.length() == n_predicted))
00535     {
00536         // No predictor part provided: this means this is the same as before
00537         // (or that there is none at all).
00538         predicted_part << input;
00539     } else {
00540         PLASSERT( input.length() == n_predictor + n_predicted );
00541         predicted_part << input.subVec(n_predictor, n_predicted);
00542         setPredictor(input);
00543     }
00544 }
00545 
00547 // forget //
00549 void PDistribution::forget() {
00550     stage = 0;
00551     n_predictor = -1;
00552     n_predicted = -1;
00553     resetGenerator(seed_);
00554 }
00555 
00557 // subclass stuff //
00559 
00560 real PDistribution::log_density(const Vec& y) const
00561 { PLERROR("density not implemented for this PDistribution"); return 0; }
00562 
00563 real PDistribution::density(const Vec& y) const
00564 { return exp(log_density(y)); }
00565 
00566 real PDistribution::survival_fn(const Vec& y) const
00567 { PLERROR("survival_fn not implemented for this PDistribution"); return 0; }
00568 
00569 real PDistribution::cdf(const Vec& y) const
00570 { PLERROR("cdf not implemented for this PDistribution"); return 0; }
00571 
00572 void PDistribution::expectation(Vec& mu) const
00573 { PLERROR("expectation not implemented for this PDistribution"); }
00574 
00575 void PDistribution::missingExpectation(const Vec& input, Vec& mu)
00576 { PLERROR("missingExpectation not implemented for this PDistribution"); }
00577 
00578 void PDistribution::variance(Mat& covar) const
00579 { PLERROR("variance not implemented for this PDistribution"); }
00580 
00581 void PDistribution::generate(Vec& y) const
00582 { PLERROR("generate not implemented for this PDistribution"); }
00583 
00584 void PDistribution::generateJoint(Vec& xy)
00585 {
00586     // get old sizes
00587     int old_n_predictor = n_predictor;
00588     int old_n_predicted = n_predicted;
00589 
00590     // set all inputs as predicted to generate a joint sample
00591     setPredictorPredictedSizes(0, -1);
00592     generate( xy );
00593 
00594     // restore old sizes
00595     setPredictorPredictedSizes(old_n_predictor, old_n_predicted);
00596 }
00597 
00598 void PDistribution::generateJoint(Vec& x, Vec& y)
00599 {
00600     Vec joint_sample;
00601     generateJoint( joint_sample );
00602     x = joint_sample.subVec(0, n_predictor);
00603     y = joint_sample.subVec(n_predictor, n_predicted);
00604 }
00605 
00606 void PDistribution::generatePredictor(Vec& x)
00607 {
00608     Vec y;
00609     generateJoint(x, y);
00610 }
00611 
00612 void PDistribution::generatePredicted(Vec& y)
00613 {
00614     Vec x;
00615     generateJoint(x, y);
00616 }
00617 
00618 void PDistribution::generatePredictorGivenPredicted(Vec& x, const Vec& y)
00619 { PLERROR("generatePredictorGivenPredicted not implemented for this\n"
00620           "PDistribution\n"); }
00621 
00622 void PDistribution::train()
00623 { PLERROR("The train() method is not implemented for this PDistribution"); }
00624 
00626 // unknownOutput //
00628 void PDistribution::unknownOutput(char def, const Vec& input, Vec& output,
00629                                   int& k) const
00630 {
00631     // Default is to throw an error.
00632     // TODO Can we find a better way to do this?
00633     PLERROR("In PDistribution::unknownOutput - Unrecognized outputs_def "
00634             "character: '%c'", def);
00635 }
00636 
00637 } // end of namespace PLearn
00638 
00639 
00640 /*
00641   Local Variables:
00642   mode:c++
00643   c-basic-offset:4
00644   c-file-style:"stroustrup"
00645   c-file-offsets:((innamespace . 0)(inline-open . 0))
00646   indent-tabs-mode:nil
00647   fill-column:79
00648   End:
00649 */
00650 // 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