/*
 * $Id: example4.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 <string>
#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;

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>&);

// generic, but simple, factorial
template <typename IntegralType>
IntegralType
factorial (IntegralType i)
{
  if (i > 0)
    return i * factorial(i-1);
  return IntegralType(1);
}

// generic, but simple, Poisson distribution
template <typename Probability, typename IntegralType>
Probability
poisson (IntegralType i)
{
  typedef typename Probability::value_type value_type;
  const value_type lambda (2);
  return Probability(exp(-lambda) * pow(lambda, i) / factorial(i),
		     linear_domain());
}

// calculate the likelihood of a set of observations
template <typename Likelihood, typename Observations>
Likelihood
analyze (const Observations& obs)
{
  Likelihood l;
  // product of likelihoods across a series of independent observations
  for (typename Observations::const_iterator i = obs.begin(); i != obs.end(); ++i)
    l *= poisson<probability>(*i);
  return l;
}

// output a label and a likelihood
template <typename Likelihood>
void
output (const std::string& label, const Likelihood& l)
{
  std::cout << label << "  " << l << std::endl;
}

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

  output("underflow:   ",analyze<likelihood,observations>(obs));
  output("no underflow:",analyze<log_likelihood,observations>(obs));

  return 0;
}

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);
}