JLegendreEstimator.hh

Go to the documentation of this file.

#ifndef __JFIT__JLEGENDREESTIMATOR__
#define __JFIT__JLEGENDREESTIMATOR__
 
#include "JMath/JLegendre.hh"
#include "JMath/JZero.hh"
#include "JMath/JMatrixNS.hh"
#include "JFit/JEstimator.hh"
 
 
/**
 * \file
 * Linear fit of JMATH::JLegendre.
 * \author mdejong
 */
namespace JFIT {}
namespace JPP { using namespace JFIT; }
 
namespace JFIT {
 
  using JMATH::JLegendre;
 
  /**
   * Linear fit of Legendre polynomial.
   */
  template<class JOrdinate_t, size_t N>
  class JEstimator< JLegendre<JOrdinate_t, N> > :
    public JLegendre<JOrdinate_t, N>
  {
  public:
    typedef JOrdinate_t                    ordinate_type;
 
 
    /**
     * Constructor.
     *
     * The template argument <tt>T</tt> refers to an iterator of a data structure which should have the following member methods:
     *   - double        %getX();       //
     *   - ordinate_type %getY();       //
     * where <tt>ordinate_type</tt> refers to the first template atgument of class JLegendre.
     *
     * \param  __begin          begin of data
     * \param  __end            end   of data
     */
    template<class T>
    JEstimator(T __begin, T __end)
    {
      using namespace std;
      using namespace JPP;
 
      if (distance(__begin, __end) >= NUMBER_OF_PARAMETERS) {
 
        vector<ordinate_type> Y(NUMBER_OF_PARAMETERS, getZero<ordinate_type>());
        vector<double>        h(NUMBER_OF_PARAMETERS, getZero<double>());
 
        V.resize(NUMBER_OF_PARAMETERS);
        V.reset();
 
        T i = __begin;
 
        this->xmin = i->getX();  advance(i, distance(__begin, __end) - 1);
        this->xmax = i->getX();
 
        for (i = __begin; i != __end; ++i) {
 
          const double z = this->getX(i->getX());
 
          for (size_t n = 0; n <= N; ++n) {
            h[n] = legendre(n, z);
          }
 
          for (size_t row = 0; row <= N; ++row) {
 
            V(row, row) += h[row] * h[row];
            Y[row]      += h[row] * i->getY();
 
            for (size_t col = 0; col != row; ++col) {
              V(row, col) += h[row] * h[col];
            }
          }
        }
 
        for (size_t row = 0; row <= N; ++row) {
          for (size_t col = 0; col != row; ++col) {
            V(col, row) = V(row, col);
          }
        }
 
        V.invert();
 
        for (size_t row = 0; row <= N; ++row) {
 
          (*this)[row] = zero;
 
          for (size_t col = 0; col <= N; ++col) {
            (*this)[row] += V(row, col) * Y[col];
          }
        }
 
      } else {
        throw JValueOutOfRange("JEstimator<JLegendre>::JEstimator(): Not enough data points.");
      }
    }
 
 
    static const int NUMBER_OF_PARAMETERS  =  N + 1;  //!< number of parameters of fit
    JMATH::JMatrixNS V;                               //!< co-variance matrix of fit parameters
  };
}
 
#endif