source: branches/0.4-stable/yat/classifier/NBC.cc @ 1392

// $Id: NBC.cc 1392 2008-07-28 19:35:30Z peter $

/*
  Copyright (C) 2006, 2007 Jari Häkkinen, Peter Johansson, Markus Ringnér
  Copyright (C) 2008 Peter Johansson, Markus Ringnér

  This file is part of the yat library, http://dev.thep.lu.se/yat

  The yat library is free software; you can redistribute it and/or
  modify it under the terms of the GNU General Public License as
  published by the Free Software Foundation; either version 2 of the
  License, or (at your option) any later version.

  The yat library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "NBC.h"
#include "MatrixLookup.h"
#include "MatrixLookupWeighted.h"
#include "Target.h"
#include "yat/statistics/Averager.h"
#include "yat/statistics/AveragerWeighted.h"
#include "yat/utility/Matrix.h"

#include <cassert>
#include <cmath>
#include <stdexcept>
#include <vector>

namespace theplu {
namespace yat {
namespace classifier {

  NBC::NBC() 
    : SupervisedClassifier()
  {
  }


  NBC::~NBC()   
  {
  }


  NBC* NBC::make_classifier() const 
  {     
    return new NBC();
  }


  void NBC::train(const MatrixLookup& data, const Target& target)
  {   
    sigma2_.resize(data.rows(), target.nof_classes());
    centroids_.resize(data.rows(), target.nof_classes());
    
    for(size_t i=0; i<data.rows(); ++i) {
      std::vector<statistics::Averager> aver(target.nof_classes());
      for(size_t j=0; j<data.columns(); ++j) 
        aver[target(j)].add(data(i,j));
      
      assert(centroids_.columns()==target.nof_classes());
      for (size_t j=0; j<target.nof_classes(); ++j){
        assert(i<centroids_.rows());
        assert(j<centroids_.columns());
        assert(i<sigma2_.rows());
        assert(j<sigma2_.columns());
        if (aver[j].n()>1){
          sigma2_(i,j) = aver[j].variance();
          centroids_(i,j) = aver[j].mean();
        }
        else {
            sigma2_(i,j) = std::numeric_limits<double>::quiet_NaN();
            centroids_(i,j) = std::numeric_limits<double>::quiet_NaN();
        }
      }
    }
  }   


  void NBC::train(const MatrixLookupWeighted& data, const Target& target)
  {   
    sigma2_.resize(data.rows(), target.nof_classes());
    centroids_.resize(data.rows(), target.nof_classes());

    for(size_t i=0; i<data.rows(); ++i) {
      std::vector<statistics::AveragerWeighted> aver(target.nof_classes());
      for(size_t j=0; j<data.columns(); ++j) 
        aver[target(j)].add(data.data(i,j), data.weight(i,j));
      
      assert(centroids_.columns()==target.nof_classes());
      for (size_t j=0; j<target.nof_classes(); ++j) {
        assert(i<centroids_.rows());
        assert(j<centroids_.columns());
        assert(i<sigma2_.rows());
        assert(j<sigma2_.columns());
        if (aver[j].n()>1){
          sigma2_(i,j) = aver[j].variance();
          centroids_(i,j) = aver[j].mean();
        }
        else {
          sigma2_(i,j) = std::numeric_limits<double>::quiet_NaN();
          centroids_(i,j) = std::numeric_limits<double>::quiet_NaN();
        }
      }
    }
  }


  void NBC::predict(const MatrixLookup& ml,                     
                    utility::Matrix& prediction) const
  {   
    assert(ml.rows()==sigma2_.rows());
    assert(ml.rows()==centroids_.rows());
    // each row in prediction corresponds to a sample label (class)
    prediction.resize(centroids_.columns(), ml.columns(), 0);

    // first calculate -lnP = sum sigma_i + (x_i-m_i)^2/2sigma_i^2
    for (size_t label=0; label<centroids_.columns(); ++label) {
      double sum_log_sigma = sum_logsigma(label);
      for (size_t sample=0; sample<prediction.rows(); ++sample) {
        prediction(label,sample) = sum_log_sigma;
        for (size_t i=0; i<ml.rows(); ++i) 
          prediction(label, sample) += 
            std::pow(ml(i, label)-centroids_(i, label),2)/
            sigma2_(i, label);
      }
    }
    standardize_lnP(prediction);
  }

      
  void NBC::predict(const MatrixLookupWeighted& mlw,                    
                    utility::Matrix& prediction) const
  {   
    assert(mlw.rows()==sigma2_.rows());
    assert(mlw.rows()==centroids_.rows());
    
    // each row in prediction corresponds to a sample label (class)
    prediction.resize(centroids_.columns(), mlw.columns(), 0);

    // first calculate -lnP = sum (sigma_i) + 
    // N sum w_i(x_i-m_i)^2/2sigma_i^2 / sum w_i
    for (size_t label=0; label<centroids_.columns(); ++label) {
      double sum_log_sigma = sum_logsigma(label);
      for (size_t sample=0; sample<prediction.rows(); ++sample) {
        statistics::AveragerWeighted aw;
        for (size_t i=0; i<mlw.rows(); ++i) 
          aw.add(std::pow(mlw.data(i, label)-centroids_(i, label),2)/
                 sigma2_(i, label), mlw.weight(i, label));
        prediction(label,sample) = sum_log_sigma + mlw.rows()*aw.mean()/2;
      }
    }
    standardize_lnP(prediction);
  }

  void NBC::standardize_lnP(utility::Matrix& prediction) const
  {
    /// -lnP might be a large number, in order to avoid out of bound
    /// problems when calculating P = exp(- -lnP), we centralize matrix
    /// by adding a constant.
    // lookup of prediction with zero weights for NaNs
    MatrixLookupWeighted mlw(prediction);
    statistics::AveragerWeighted a;
    add(a, mlw.begin(), mlw.end());
    prediction -= a.mean();
    
    // exponentiate
    for (size_t i=0; i<prediction.rows(); ++i)
      for (size_t j=0; j<prediction.columns(); ++j)
        prediction(i,j) = std::exp(prediction(i,j));
    
    // normalize each row (label) to sum up to unity (probability)
    for (size_t i=0; i<prediction.rows(); ++i){
      // calculate sum of row ignoring NaNs
      statistics::AveragerWeighted a;
      add(a, mlw.begin_row(i), mlw.end_row(i));
      prediction.row_view(i) *= 1.0/a.sum_wx();
    }
  }


  double NBC::sum_logsigma(size_t label) const
  {
    double sum_log_sigma=0;
    assert(label<sigma2_.columns());
    for (size_t i=0; i<sigma2_.rows(); ++i) {
      sum_log_sigma += std::log(sigma2_(i, label));
    }
    return sum_log_sigma / 2; // taking sum of log(sigma) not sigma2
  }

}}} // of namespace classifier, yat, and theplu