/*
 * $Id: example1.cpp,v 1.2 2007/06/27 16:38:24 brook Exp $
 */

/*
 * Copyright (c) 2007 Brook Milligan.
 * Distributed under the Boost Software License, Version 1.0. (See
 * accompanying file LICENSE_1_0.txt or copy at
 * http://www.boost.org/LICENSE_1_0.txt)
 */

#if defined _MSC_VER
#  pragma warning(disable: 4100) // 'b' : unreferenced formal parameter
#endif

#include <cmath>
#include <iostream>
#include <vector>
#include <boost/likelihood.hpp>
#include <boost/probability.hpp>

typedef boost::probabilities::linear_domain linear_domain;
typedef boost::probabilities::log_domain log_domain;

typedef boost::likelihood likelihood;
typedef boost::log_likelihood log_likelihood;
typedef boost::probability probability;

typedef std::vector<int> observations;

static probability poisson (int i); // probability model

template <typename Value, typename Validator>
static std::ostream& operator << (std::ostream&, const boost::probabilities::
				  likelihood<linear_domain,Value,Validator>&);
template <typename Value, typename Validator>
static std::ostream& operator << (std::ostream&, const boost::probabilities::
				  likelihood<log_domain,Value,Validator>&);

int
main ()
{
  // a series of observations
  observations obs;
  for (int i = 0; i <= 10; ++i)
    obs.push_back(i % 10);   // limit observations to the range [0, 9]

  log_likelihood l;
  // product of likelihoods across a series of independent observations
  for (observations::const_iterator i = obs.begin(); i != obs.end(); ++i)
    l *= poisson(*i);

  std::cout << boost::probabilities::domain_cast<linear_domain>(l)
	    << std::endl;
  std::cout << l << std::endl;

  return 0;
}

int
factorial (int i)
{
  if (i > 0)
    return i * factorial(i-1);
  return 1;
}

probability
poisson (int i)
{
  const double lambda (2);
  return probability(exp(-lambda) * pow(lambda, i) / factorial(i));
}

template <typename Value, typename Validator>
std::ostream&
operator << (std::ostream& s, const boost::probabilities::
	     likelihood<linear_domain,Value,Validator>& l)
{
  return s << "l=" << boost::probabilities::value_cast<linear_domain>(l);
}

template <typename Value, typename Validator>
std::ostream&
operator << (std::ostream& s, const boost::probabilities::
	     likelihood<log_domain,Value,Validator>& l)
{
  return s << "ln(l)=" << boost::probabilities::value_cast<log_domain>(l);
}