source: trunk/yat/classifier/Perceptron.cc @ 3709

// $Id: Perceptron.cc 3709 2017-11-08 22:49:06Z peter $

/*
  Copyright (C) 2017 Peter Johansson

  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 3 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 yat. If not, see <http://www.gnu.org/licenses/>.
*/

#include <config.h>

#include "Perceptron.h"

#include "Target.h"

#include "yat/utility/DiagonalMatrix.h"
#include "yat/utility/Matrix.h"
#include "yat/utility/Vector.h"

#include <gsl/gsl_cdf.h>

#include <cassert>
#include <cmath>
#include <cmath>

namespace theplu {
namespace yat {
namespace classifier {

  const utility::Matrix& Perceptron::covariance(void) const
  {
    return covariance_;
  }


  double Perceptron::margin(size_t i, double alpha) const
  {
    return gsl_cdf_ugaussian_Qinv(alpha/2) * std::sqrt(covariance_(i, i));
  }


  double Perceptron::oddsratio(size_t i) const
  {
    return std::exp(weight_(i));
  }


  double Perceptron::oddsratio_lower_CI(size_t i, double alpha) const
  {
    return std::exp(weight_(i) - margin(i, alpha));
  }


  double Perceptron::oddsratio_upper_CI(size_t i, double alpha) const
  {
    return std::exp(weight_(i) + margin(i, alpha));
  }


  double Perceptron::p_value(size_t i) const
  {
    double z = weight_(i) / std::sqrt(covariance_(i, i));
    return 2*gsl_cdf_ugaussian_Q(std::abs(z));
  }


  double Perceptron::predict(const utility::VectorBase& x) const
  {
    assert(x.size() == weight_.size());
    const double f = weight_ * x;
    return 1.0 / (1 + std::exp(-f));
  }


  void Perceptron::train(const utility::Matrix& X, const Target& target)
  {
    size_t n = X.rows();
    size_t p = X.columns();

    assert(target.size() == n);
    weight_.resize(p);
    covariance_.resize(p, p);

    // weight vector is updated as
    // w = (X'SX)^-1 X' (SXw + y - mu)
    // X is n x p
    // mu is vector of (trained) expected values (see predict(1))
    utility::Vector mu(n);
    // S is diagonal n x n with S_ii = mu_i (1 - mu_i)
    utility::DiagonalMatrix S(n, n);
    // y is binary vector
    utility::Vector y(n);
    for (size_t i=0; i<n; ++i)
      if (target.binary(i))
        y(i) = 1.0;

    size_t max_epochs = 100;
    double sum_squared = 1.0; // some (relatively) large number
    for (size_t epoch=0; sum_squared > 1e-8 && epoch < max_epochs; ++epoch) {
      for (size_t i=0; i<mu.size(); ++i) {
        mu(i) = predict(X.row_const_view(i));
        S(i) = mu(i) * (1.0 - mu(i));
      }

      // w = (X'SX)^-1 X' (SXw + y - mu)
      assert(X.rows() == S.rows());
      assert(S.columns() == X.rows());
      utility::inverse_svd(transpose(X)*S*X, covariance_);

      assert(y.size() == mu.size());
      utility::Vector delta = covariance_ * (transpose(X) * (y - mu));
      weight_ += delta;
      sum_squared = delta * delta;
    }
  }


  const utility::Vector& Perceptron::weight(void) const
  {
    return weight_;
  }


}}}