source: trunk/test/poisson.cc

// $Id: poisson.cc 3662 2017-07-19 01:53:55Z 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 "Suite.h"

#include "yat/regression/Poisson.h"
#include "yat/random/random.h"
#include "yat/statistics/Averager.h"
#include "yat/utility/Matrix.h"
#include "yat/utility/Vector.h"

#include <vector>

using namespace theplu::yat;

void analyse(const utility::Vector& b, const utility::Matrix& X,
             std::vector<statistics::Averager>& stats,
             std::vector<statistics::Averager>& stats2);

void generate_X(const utility::Vector& b, utility::Matrix& X);

int main(int argc, char* argv[])
{
  test::Suite suite(argc, argv);

  regression::Poisson model;
  size_t n = 1000;
  size_t p = 4;
  utility::Vector y(n);
  utility::Matrix X(n, p);
  model.fit(X, y);
  if (model.fit_parameters().size() != p+1) {
    suite.add(false);
    suite.err() << "error: size of fit_parameters: "
                << model.fit_parameters().size()
                << " expected " << p+1 << "\n";
  }
  model.predict(X.row_const_view(0));

  utility::Vector b(4);
  b(0) = 0.5;
  b(1) = 2.0;
  b(2) = 0.75;
  b(3) = -1.25;
  generate_X(b, X);
  std::vector<statistics::Averager> stats(b.size());
  std::vector<statistics::Averager> stats2(b.size());
  for (size_t i=0; i<100; ++i)
    analyse(b, X, stats, stats2);

  for (size_t i=0; i<b.size(); ++i) {
    suite.out() << i << " " << b(i) << " " << stats[i].mean() << " "
                << stats[i].standard_error() << "\n";
    if (stats[i].standard_error() == 0.0) {
      suite.add(false);
      suite.err() << "error: standard error is 0.0\n";
    }
    else if (std::abs(stats[i].mean() - b(i)) > 5*stats[i].standard_error()) {
      suite.add(false);
      suite.err() << "error: average for param " << i << ": "
                  << stats[i].mean() << "\n";
    }
  }

  suite.out() << "\ntest covariance\n";
  suite.out() << "i\t<inferred>\tobserved\tdelta\trelative error\n";
  for (size_t i=0; i<b.size(); ++i) {
    // comparing uncertainty observed (stats2) and uncertainty
    // estimated (stats)

    double delta = stats2[i].mean() - stats[i].std();
    double relative_error = delta / stats[i].std();

    suite.out() << i << "\t" << stats2[i].mean() << "\t"
                << stats[i].std() << "\t"
                << delta << "\t"
                << relative_error << "\n";

    // covariance matrix is only an approximation, so only check that
    // estimation is roughly correct (within 20%)
    if (std::abs(relative_error) > 0.2) {
      suite.add(false);
      suite.err() << "error: " << relative_error << " too large\n";
    }
  }

  return suite.return_value();
}


void generate_X(const utility::Vector& b, utility::Matrix& X)
{
  size_t n = 5000;
  X.resize(n, b.size()-1);
  random::Gaussian gauss;
  for (utility::Matrix::iterator it=X.begin(); it!=X.end(); ++it)
    *it = gauss();
}


void generate_y(const utility::Vector& b,
                const utility::Matrix& X, utility::Vector& y)
{
  size_t n = X.rows();
  y.resize(n);
  random::Poisson poisson;

  for (size_t i=0; i<n; ++i) {
    double lnmu = b(0);
    for (size_t j=1; j<b.size(); ++j)
      lnmu += X(i, j-1) * b(j);

    double mu = exp(lnmu);
    y(i) = poisson(mu);
  }
}


void analyse(const utility::Vector& b, const utility::Matrix& X,
             std::vector<statistics::Averager>& stats,
             std::vector<statistics::Averager>& stats2)
{
  utility::Vector y;
  generate_y(b, X, y);
  theplu::yat::regression::Poisson poisson;
  poisson.fit(X, y);
  for (size_t i=0; i<poisson.fit_parameters().size(); ++i) {
    stats[i].add(poisson.fit_parameters()(i));
    stats2[i].add(std::sqrt(poisson.covariance()(i, i)));
  }
}