Merge branch 'issue#593' of https://github.com/ShihaoJing/mlpack into ShihaoJing-issue#593

src/mlpack/methods/naive_bayes/naive_bayes_classifier.hpp

@@ -1,6 +1,7 @@
 /**
  * @file naive_bayes_classifier.hpp
  * @author Parikshit Ram (pram@cc.gatech.edu)
+ * @author Shihao Jing (shihao.jing810@gmail.com)
  *
  * A Naive Bayes Classifier which parametrically estimates the distribution of
  * the features.  It is assumed that the features have been sampled from a
@@ -110,20 +111,67 @@ class NaiveBayesClassifier
   void Train(const VecType& point, const size_t label);
 
   /**
-   * Given a bunch of data points, this function evaluates the class of each of
-   * those data points, and puts it in the vector 'results'.
+   * Classify the given point, using the training GaussianNB model. The predicted label is
+   * returned.
+   *
+   * @param point Point to classify.
+   */
+  template<typename VecType>
+  size_t Classify(const VecType& point) const;
+
+  /**
+   * Classify the given point using the training GaussianNB model 
+   * and also return estimates of the probability for
+   * each class in the given vector.
+   *
+   * @param point Point to classify.
+   * @param prediction This will be set to the predicted class of the point.
+   * @param probabilities This will be filled with class probabilities for the
+   *      point.
+   */
+  template<typename VecType>
+  void Classify(const VecType& point,
+                size_t& prediction,
+                arma::vec& probabilities) const;
+
+   /**
+   * Classify the given points using the training GaussianNB model.
+   * The predicted labels for each point are stored in the given vector.
    *
    * @code
    * arma::mat test_data; // each column is a test point
    * arma::Row<size_t> results;
    * ...
-   * nbc.Classify(test_data, &results);
+   * nbc.Classify(test_data, results);
    * @endcode
    *
    * @param data List of data points.
-   * @param results Vector that class predictions will be placed into.
+   * @param predictions that class predictions will be placed into.
    */
-  void Classify(const MatType& data, arma::Row<size_t>& results);
+  void Classify(const MatType& data, 
+                arma::Row<size_t>& predictions) const;
+
+  /**
+   * Classify the given points using the training GaussianNB model 
+   * and also return estimates of the probabilities for each class in the given matrix.  
+   * The predicted labels for each point are stored in the given vector.
+   *
+   * @code
+   * arma::mat test_data; // each column is a test point
+   * arma::Row<size_t> results;
+   * arma::mat resultsProbs;
+   * ...
+   * nbc.Classify(test_data, results, resultsProbs);
+   * @endcode
+   *
+   * @param data Set of points to classify.
+   * @param predictions This will be filled with predictions for each point.
+   * @param probabilities This will be filled with class probabilities for each
+   *      point. Each row represents a point.
+   */
+  void Classify(const MatType& data,
+                arma::Row<size_t>& predictions,
+                arma::mat& probabilities) const;
 
   //! Get the sample means for each class.
   const MatType& Means() const { return means; }
@@ -145,6 +193,7 @@ class NaiveBayesClassifier
   void Serialize(Archive& ar, const unsigned int /* version */);
 
  private:
+
   //! Sample mean for each class.
   MatType means;
   //! Sample variances for each class.
@@ -153,6 +202,23 @@ class NaiveBayesClassifier
   arma::vec probabilities;
   //! Number of training points seen so far.
   size_t trainingPoints;
+
+  /**
+  * Compute the unnormalized posterior log probability of given point.
+  *
+  * @param point Data point to compute posterior log probability. 
+  */
+  template<typename VecType>
+  arma::vec JointLogLikelihood(const VecType& point) const;
+
+  /**
+  * Compute the unnormalized posterior log probability of given points.
+  * Results are returned as arma::mat, and each ling represents a points,
+  * each column represents posterior log probability of a class.
+  *
+  * @param @param data Set of points to compute posterior log probability. 
+  */
+  arma::mat JointLogLikelihood(const MatType& point) const;
 };
 
 } // namespace naive_bayes

src/mlpack/methods/naive_bayes/naive_bayes_classifier_impl.hpp

@@ -2,6 +2,7 @@
  * @file naive_bayes_classifier_impl.hpp
  * @author Parikshit Ram (pram@cc.gatech.edu)
  * @author Vahab Akbarzadeh (v.akbarzadeh@gmail.com)
+ * @author Shihao Jing (shihao.jing810@gmail.com)
  *
  * A Naive Bayes Classifier which parametrically estimates the distribution of
  * the features.  This classifier makes its predictions based on the assumption
@@ -162,24 +163,18 @@ void NaiveBayesClassifier<MatType>::Train(const VecType& point,
 }
 
 template<typename MatType>
-void NaiveBayesClassifier<MatType>::Classify(const MatType& data,
-                                             arma::Row<size_t>& results)
+template<typename VecType>
+arma::vec NaiveBayesClassifier<MatType>::JointLogLikelihood(
+                                      const VecType& point) const
 {
   // Check that the number of features in the test data is same as in the
   // training data.
-  Log::Assert(data.n_rows == means.n_rows);
+  Log::Assert(point.n_rows == means.n_rows);
 
-  arma::vec probs = arma::log(probabilities);
+  arma::vec jll = arma::log(probabilities);
   arma::mat invVar = 1.0 / variances;
 
-  arma::mat testProbs = arma::repmat(probs.t(), data.n_cols, 1);
-
-  results.set_size(data.n_cols); // No need to fill with anything yet.
-
-  Log::Info << "Running Naive Bayes classifier on " << data.n_cols
-      << " data points with " << data.n_rows << " features each." << std::endl;
-
-  // Calculate the joint probability for each of the data points for each of the
+  // Calculate the joint log likelihood of point for each of the
   // means.n_cols.
 
   // Loop over every class.
@@ -187,30 +182,101 @@ void NaiveBayesClassifier<MatType>::Classify(const MatType& data,
   {
     // This is an adaptation of gmm::phi() for the case where the covariance is
     // a diagonal matrix.
-    arma::mat diffs = data - arma::repmat(means.col(i), 1, data.n_cols);
-    arma::mat rhs = -0.5 * arma::diagmat(invVar.col(i)) * diffs;
-    arma::vec exponents(diffs.n_cols);
-    for (size_t j = 0; j < diffs.n_cols; ++j) // log(exp(value)) == value
-      exponents(j) = arma::accu(diffs.col(j) % rhs.unsafe_col(j));
+    arma::vec diffs = point - means.col(i);
+    arma::vec rhs = -0.5 * arma::diagmat(invVar.col(i)) * diffs;
+    double exponent = arma::accu(diffs % rhs); // log(exp(value)) == value
 
-    // Calculate probability as sum of logarithm to decrease floating point
-    // errors.
-    testProbs.col(i) += (data.n_rows / -2.0 * log(2 * M_PI) - 0.5 *
-        log(arma::det(arma::diagmat(variances.col(i)))) + exponents);
+    // Calculate oint log likelihood as sum of logarithm 
+    // to decrease floating point errors.
+    jll(i) += (point.n_rows / -2.0 * log(2 * M_PI) - 0.5 *
+        log(arma::det(arma::diagmat(variances.col(i)))) + exponent);
   }
 
-  // Now calculate the label.
+  return jll;
+}
+
+template<typename MatType>
+arma::mat NaiveBayesClassifier<MatType>::JointLogLikelihood(
+                                      const MatType& data) const
+{
+  arma::mat jll(data.n_cols, means.n_cols);
   for (size_t i = 0; i < data.n_cols; ++i)
   {
-    // Find the index of the class with maximum probability for this point.
-    arma::uword maxIndex = 0;
-    arma::vec pointProbs = testProbs.row(i).t();
-    pointProbs.max(maxIndex);
-
-    results[i] = maxIndex;
+    jll.row(i) = JointLogLikelihood(data.col(i)).t();
   }
+  return jll;
+}        
 
-  return;
+template<typename MatType>
+template<typename VecType>
+size_t NaiveBayesClassifier<MatType>::Classify(const VecType& point) const
+{
+  // find the label(class) with max joint log likelihood.
+  arma::vec jll = JointLogLikelihood(point);
+  arma::uword maxIndex = 0;
+  jll.max(maxIndex);
+  return maxIndex;
+}
+
+template<typename MatType>
+template<typename VecType>
+void NaiveBayesClassifier<MatType>::Classify(const VecType& point,
+                                             size_t& prediction,
+                                             arma::vec& probabilities) const
+{
+  // log(Prob(Y|X)) = Log(X|Y) + Log(Y) - Log(X);
+  // JointLogLikelihood = Log(X|Y) + Log(Y)
+  arma::vec jll = JointLogLikelihood(point);
+  double logProbX = log(arma::accu(exp(jll))); // Log(X)
+  jll -= logProbX;
+
+  arma::uword maxIndex = 0;
+  jll.max(maxIndex);
+  prediction = maxIndex;
+  probabilities = exp(jll); // log(exp(value)) == value
+}
+
+template<typename MatType>
+void NaiveBayesClassifier<MatType>::Classify(
+                                 const MatType& data,
+                                 arma::Row<size_t>& predictions) const
+{
+  predictions.set_size(data.n_cols);
+
+  arma::mat jll = JointLogLikelihood(data);
+
+  for (size_t i = 0; i < data.n_cols; ++i)
+  {
+    arma::uword maxIndex = 0;
+    arma::vec pointProbs = jll.row(i).t();
+    pointProbs.max(maxIndex);
+    predictions[i] = maxIndex;
+  }
+}
+
+template<typename MatType>
+void NaiveBayesClassifier<MatType>::Classify(const MatType& data,
+                arma::Row<size_t>& predictions,
+                arma::mat& predictionProbs) const
+{
+  predictions.set_size(data.n_cols);
+
+  arma::mat jll = JointLogLikelihood(data);
+  arma::vec logProbX(data.n_cols); // log(Prob(X))
+  for (size_t j = 0; j < data.n_cols; ++j)
+    logProbX(j) = log(arma::accu(exp(jll.row(j))));
+
+  jll -= arma::repmat(logProbX, 1, data.n_cols);
+
+  predictionProbs = exp(jll);
+
+  for (size_t i = 0; i < data.n_cols; ++i)
+  {
+    arma::uword maxIndex = 0;
+    arma::vec pointProbs = jll.row(i).t();
+    pointProbs.max(maxIndex);
+    predictions[i] = maxIndex;
+  }
 }
 
 template<typename MatType>

src/mlpack/methods/naive_bayes/nbc_main.cpp

@@ -77,6 +77,8 @@ PARAM_FLAG("incremental_variance", "The variance of each class will be "
 PARAM_MATRIX_IN("test", "A matrix containing the test set.", "T");
 PARAM_UROW_OUT("output", "The matrix in which the predicted labels for the"
     " test set will be written.", "o");
+PARAM_MATRIX_OUT("output_probs", "The matrix in which the predicted probability of labels for the"
+    " test set will be written.", "p");
 
 int main(int argc, char* argv[])
 {
@@ -98,8 +100,8 @@ int main(int argc, char* argv[])
     Log::Warn << "--incremental_variance (-I) ignored because --training_file "
         << "(-t) is not specified." << endl;
 
-  if (!CLI::HasParam("output") && !CLI::HasParam("output_model"))
-    Log::Warn << "Neither --output_file (-o) nor --output_model_file (-M) "
+  if (!CLI::HasParam("output") && !CLI::HasParam("output_model") && !CLI::HasParam("output_probs"))
+    Log::Warn << "Neither --output_file (-o), nor --output_model_file (-M), nor --output_proba_file (-p)"
         << "specified; no output will be saved!" << endl;
 
   if (CLI::HasParam("output") && !CLI::HasParam("test"))
@@ -110,6 +112,10 @@ int main(int argc, char* argv[])
     Log::Warn << "--test_file (-T) specified, but classification results will "
         << "not be saved because --output_file (-o) is not specified." << endl;
 
+  if (!CLI::HasParam("output_probs") && CLI::HasParam("test"))
+    Log::Warn << "--test_file (-T) specified, but predicted probability of labels will "
+        << "not be saved because --output_probs_file (-p) is not specified." << endl;
+
   // Either we have to train a model, or load a model.
   NBCModel model;
   if (CLI::HasParam("training"))
@@ -160,20 +166,24 @@ int main(int argc, char* argv[])
           << ")!" << std::endl;
 
     // Time the running of the Naive Bayes Classifier.
-    Row<size_t> results;
+    Row<size_t> predictions;
+    mat probabilities;
     Timer::Start("nbc_testing");
-    model.nbc.Classify(testingData, results);
+    model.nbc.Classify(testingData, predictions, probabilities);
     Timer::Stop("nbc_testing");
 
     if (CLI::HasParam("output"))
     {
       // Un-normalize labels to prepare output.
       Row<size_t> rawResults;
-      data::RevertLabels(results, model.mappings, rawResults);
+      data::RevertLabels(predictions, model.mappings, rawResults);
 
       // Output results.
       CLI::GetParam<Row<size_t>>("output") = std::move(rawResults);
     }
+
+    if (CLI::HasParam("output_probs"))
+      CLI::GetParam<mat>("output_probs") = probabilities.t();
   }
 
   if (CLI::HasParam("output_model"))

src/mlpack/tests/data/testResProbs.csv

@@ -0,0 +1,7 @@
+0.999999999998,1.86113163315e-12
+0.999999999996,3.82159650359e-12
+0.999999999998,2.42337174628e-12
+2.96941920603e-19,1.0
+2.96941920603e-19,1.0
+1.44611813707e-19,1.0
+2.28049205613e-19,1.0

src/mlpack/tests/nbc_test.cpp

@@ -25,6 +25,7 @@ BOOST_AUTO_TEST_CASE(NaiveBayesClassifierTest)
   const char* testFilename = "testSet.csv";
   const char* trainResultFilename = "trainRes.csv";
   const char* testResultFilename = "testRes.csv";
+  const char* testResultProbsFilename = "testResProbs.csv";
   size_t classes = 2;
 
   arma::mat trainData, trainRes, calcMat;
@@ -60,16 +61,28 @@ BOOST_AUTO_TEST_CASE(NaiveBayesClassifierTest)
 
   arma::mat testData;
   arma::Mat<size_t> testRes;
+  arma::mat testResProbs;
   arma::Row<size_t> calcVec;
+  arma::mat calcProbs;
   data::Load(testFilename, testData, true);
   data::Load(testResultFilename, testRes, true);
+  data::Load(testResultProbsFilename, testResProbs, true);
 
   testData.shed_row(testData.n_rows - 1); // Remove the labels.
 
-  nbcTest.Classify(testData, calcVec);
+  nbcTest.Classify(testData, calcVec, calcProbs);
+  calcProbs = calcProbs.t();
 
   for (size_t i = 0; i < testData.n_cols; i++)
     BOOST_REQUIRE_EQUAL(testRes(i), calcVec(i));
+
+  for (size_t i = 0; i < testResProbs.n_cols; ++i)
+  {
+    for (size_t j = 0; j < testResProbs.n_rows; ++j)
+    {
+      BOOST_REQUIRE_CLOSE(testResProbs(i, j) + 0.0001, calcProbs(i, j) + 0.0001, 0.01);
+    }
+  }
 }
 
 // The same test, but this one uses the incremental algorithm to calculate
@@ -80,6 +93,7 @@ BOOST_AUTO_TEST_CASE(NaiveBayesClassifierIncrementalTest)
   const char* testFilename = "testSet.csv";
   const char* trainResultFilename = "trainRes.csv";
   const char* testResultFilename = "testRes.csv";
+  const char* testResultProbsFilename = "testResProbs.csv";
   size_t classes = 2;
 
   arma::mat trainData, trainRes, calcMat;
@@ -115,16 +129,24 @@ BOOST_AUTO_TEST_CASE(NaiveBayesClassifierIncrementalTest)
 
   arma::mat testData;
   arma::Mat<size_t> testRes;
+  arma::mat testResProba;
   arma::Row<size_t> calcVec;
+  arma::mat calcProbs;
   data::Load(testFilename, testData, true);
   data::Load(testResultFilename, testRes, true);
+  data::Load(testResultProbsFilename, testResProba, true);
 
   testData.shed_row(testData.n_rows - 1); // Remove the labels.
 
-  nbcTest.Classify(testData, calcVec);
+  nbcTest.Classify(testData, calcVec, calcProbs);
+  calcProbs = calcProbs.t();
 
   for (size_t i = 0; i < testData.n_cols; i++)
     BOOST_REQUIRE_EQUAL(testRes(i), calcVec(i));
+
+  for (size_t i = 0; i < testResProba.n_cols; ++i)
+    for (size_t j = 0; j < testResProba.n_rows; ++j)
+      BOOST_REQUIRE_CLOSE(testResProba(i, j) + .00001, calcProbs(i, j) + .00001, 0.01);
 }
 
 /**