source: trunk/yat/statistics/AveragerPairWeighted.h @ 892

#ifndef _theplu_yat_statistics_averagerpairweighted_
#define _theplu_yat_statistics_averagerpairweighted_

// $Id: AveragerPairWeighted.h 892 2007-09-25 19:25:26Z peter $

/*
  Copyright (C) 2005 Markus Ringnér, Peter Johansson
  Copyright (C) 2006 Jari Häkkinen, Markus Ringnér, Peter Johansson
  Copyright (C) 2007 Peter Johansson

  This file is part of the yat library, http://trac.thep.lu.se/trac/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 2 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 this program; if not, write to the Free Software
  Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  02111-1307, USA.
*/

#include "AveragerWeighted.h"

#include <cmath>

namespace theplu{
namespace yat{
namespace classifier{
  class DataLookup1D;
  class DataLookupWeighted1D;
}
namespace statistics{
  ///
  /// @brief Class for taking care of mean and covariance of two variables in
  /// a weighted manner.
  ///
  /// <a href="Statistics/index.html">Weighted Statistics document</a>
  ///
  /// If nothing else stated, each function fulfills the
  /// following:<br> <ul><li>Setting a weight to zero corresponds to
  /// removing the data point from the dataset.</li><li> Setting all
  /// weights to unity, the yields the same result as from
  /// corresponding function in AveragerPair.</li><li> Rescaling weights
  /// does not change the performance of the object.</li></ul>
  ///
  /// @see Averager AveragerWeighted AveragerPair
  ///
  class AveragerPairWeighted
  {
  public:

    ///
    /// @brief The default constructor
    ///
    AveragerPairWeighted(void);

    ///
    /// Adding a pair of data points with value \a x and \a y, and
    /// their weights. If either of the weights are zero the addition
    /// is ignored
    ///
    void  add(const double x, const double y, 
              const double wx, const double wy);

    ///
    /// Adding two sequences of data @a x and @a y. The data
    /// should be paired so \f$ x(i) \f$ is associated to \f$ y(i) \f$
    /// @a x will be treated as having all weights equal to unity
    ///
    void add(const classifier::DataLookup1D& x, 
             const classifier::DataLookupWeighted1D& y);

    ///
    /// Adding two sequences of data @a x and @a y. The data should be
    /// paired so \f$ x(i) \f$ is associated to \f$ y(i) \f$
    /// @a y will be treated as having all weights equal to unity
    ///
    void add(const classifier::DataLookupWeighted1D& x, 
             const classifier::DataLookup1D& y);

    ///
    /// Adding two sequences of weighted data @a x and @a y. The data
    /// should be paired so \f$ x(i) \f$ is associated to \f$ y(i) \f$
    ///
    void add(const classifier::DataLookupWeighted1D& x, 
             const classifier::DataLookupWeighted1D& y);

    ///
    /// Adding pair of data and corresponding pair of weights in arrays
    ///
    ///
    /// The requirements for the types T1, T2, T3 and T4 of the arrays
    /// are: operator[] returning an element and function /// size()
    /// returning the number of elements.
    ///
    template <typename T1,typename T2,typename T3,typename T4>
    void add_values(const T1& x, 
                    const T2& y, 
                    const T3& wx, 
                    const T4& wy);

    ///
    /// @brief Pearson correlation coefficient.
    ///
    /// @return \f$ \frac{\sum w_xw_y (x-m_x)(y-m_y)}{\sqrt{\sum
    /// w_xw_y (y-m_y)^2\sum w_xw_y (y-m_y)^2}} \f$ where m is
    /// calculated as \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
    ///
    double correlation(void) const;
  
    ///
    /// \f$ \frac{\sum w_xw_y (x-m_x)(y-m_y)}{\sum w_xw_y} \f$ where m
    /// is calculated as \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
    ///
    double covariance(void) const;

    /**
       @return \f$ \frac{\sum w_xw_y(x-y)^2}{\sum w_xw_y} \f$
    */
    double msd(void) const;

    ///
    /// @brief Reset everything to zero
    ///
    void reset(void);

    /**
       @return Sum of weighted squared deviation between x and y \f$
       \sum w_xw_y (x-y)^2 \f$

       \note This function is not invariant under rescaling of weights.
    */
    double sum_squared_deviation(void) const;

    ///
    /// @return \f$ \sum w_xw_y \f$
    ///
    double sum_w(void) const;

    ///
    /// @return \f$ \sum w_xw_yxy \f$
    ///
    double sum_xy(void) const;

    ///
    /// @return \f$ \sum w_xw_y (x-m_x)(y-m_y) \f$ where m is calculated as
    /// \f$ m_x = \frac {\sum w_xw_yx}{\sum w} \f$
    ///
    double sum_xy_centered(void) const;

    ///
    /// @note the weights are calculated as \f$ w = w_x * w_y \f$
    ///
    /// @return AveragerWeighted for x
    ///
    const AveragerWeighted& x_averager(void) const;

    ///
    /// @note the weights are calculated as \f$ w = w_x * w_y \f$
    ///
    /// @return AveragerWeighted for y
    ///
    const AveragerWeighted& y_averager(void) const;

  private:
    AveragerWeighted x_; // weighted averager with w = w_x*w_y
    AveragerWeighted y_; // weighted averager with w = w_x*w_y
    double wxy_;
    double w_;

  };

  /**
     \brief adding a ranges of values to AveragerPairWeighted \a ap
  */
  template <class Iter1, class Iter2>
  void add(AveragerPairWeighted& ap, Iter1 first1, Iter1 last1, Iter2 first2)
  {
    for ( ; first1 != last1; ++first1, ++first2)
      ap.add(first1.data(), first2.data(),first1.weight(),first2.weight());
  }




  template <typename T1,typename T2,typename T3,typename T4>
  void  AveragerPairWeighted::add_values(const T1& x, 
                                         const T2& y, 
                                         const T3& wx,
                                         const T4& wy)
  {
    for (size_t i=0; i<x.size(); i++) 
      add(x[i],y[i],wx[i],wy[i]);
  }

}}} // of namespace statistics, yat, and theplu

#endif