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