source: trunk/test/poisson.cc @ 3648

// $Id: poisson.cc 3648 2017-06-02 06:01:27Z 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,
             std::vector<statistics::Averager>& stats,
             std::vector<statistics::Averager>& stats2);

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

  regression::Poisson model;
  size_t n = 20;
  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;
  std::vector<statistics::Averager> stats(b.size());
  std::vector<statistics::Averager> stats2(b.size());
  for (size_t i=0; i<100; ++i)
    analyse(b, 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.xadd(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.xadd(false);
      suite.err() << "error: average for param " << i << ": "
                  << stats[i].mean() << "\n";
    }
  }

  return suite.return_value();
}


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

  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,
             std::vector<statistics::Averager>& stats,
             std::vector<statistics::Averager>& stats2)
{
  utility::Matrix X;
  utility::Vector y;
  generate_data(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));
  }
}