source: trunk/c++_tools/statistics/Histogram.h @ 675

#ifndef _theplu_statistics_histogram_
#define _theplu_statistics_histogram_

// $Id: Histogram.h 675 2006-10-10 12:08:45Z jari $

/*
  Copyright (C) The authors contributing to this file.

  This file is part of the yat library, http://lev.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 "yat/statistics/AveragerWeighted.h"

#include <string>
#include <vector>


namespace theplu {
namespace statistics {

  ///
  /// Histograms provide a convenient way of presenting the
  /// distribution of a set of data. A histogram consists of a set of
  /// bins which count the number of events falling into these
  /// bins. Currently only one dimensional histograms with uniformly
  /// spaced bins are supported.
  ///
  class Histogram
  {
  public:

    ///
    /// The default constructor.
    ///
    Histogram(void);

    ///
    /// The copy constructor.
    ///
    Histogram(const Histogram&);

    ///
    /// Construct a histogram object that covers \f$(xmin,xmax]\f$
    /// with the bin spacing \f$(xmax-xmin)/n\f$.
    ///
    Histogram(const double xmin, const double xmax, const size_t n);

    virtual ~Histogram(void);

    ///
    /// Update the histogram by adding \a weight to the bin whose
    /// range contains the observation \a x. No bins are updated when
    /// \a x lies outside the range of the histogram but the value is
    /// added to the overall integral of the histogram.
    ///
    /// @short Add an data point to the histogram.
    ///
    /// @return 0 if \a x lies within the range of the histogram, -1
    /// if \a x is smaller than the lower limit of the histogram, and
    /// similarly, 1 is returned if \a x is greater than or equal to
    /// the upper limit.
    ///
    int add(const double x,const double weight=1.0);

    ///
    /// Gives access to the AveragerWeighted object that keeps track of
    /// average of all events presented to the histogram.
    ///
    /// @short Average of all events presented to the histogram.
    ///
    /// @return A const reference to an AveragerWeighted object.
    ///
    inline const statistics::AveragerWeighted& averager_all(void) const
      { return sum_all_; }

    ///
    /// Gives access to the AveragerWeighted object that keeps track of
    /// average of events that fits within the histogram lower and
    /// upper limits. This function is equivalent to averager().
    ///
    /// @short Average of events fitting within histogram.
    ///
    /// @return A const reference to an AveragerWeighted object.
    ///
    inline const statistics::AveragerWeighted& averager_histogram(void) const
      { return sum_histogram_; }

    ///
    /// @return The number of bins in the histogram
    ///
    inline size_t nof_bins(void) const { return histogram_.size(); }

    ///
    ///  There are two ways to normalize the counts.
    ///
    /// If choice is true: The normalized count is the count in a
    /// bin divided by the total number of observations. In this
    /// case the relative counts are normalized to sum to unity (
    /// minus values outside histogram).  This is the intuitive case
    /// where the height of the histogram bar represents the
    /// proportion of the data in each class.
    ///
    /// If choice is false: The normalized count is the count in the
    /// class divided by the number of observations times the bin
    /// width. For this normalization, the area (or integral) under
    /// the histogram is equal to unity (minus the missing area
    /// corresponding to counts outside histogram). From a
    /// probabilistic point of view, this normalization results in a
    /// relative histogram that is most akin to the probability
    /// density function If you want to overlay a probability density
    /// on top of the histogram, use this normalization. Although this
    /// normalization is less intuitive (relative frequencies greater
    /// than 1 are quite permissible), it is the appropriate
    /// normalization if you are using the histogram to model a
    /// probability density function.
    ///
    /// @short Normalizing the histogram
    /// 
    void normalize(bool choice = true);

    ///
    /// @return The value in the middle of bin \a k.
    ///
    /// @note No check is done that \a k is within the size of the
    /// histogram.
    ///
    inline double observation_value(const size_t k) const
      { return xmin_+spacing()*(k+0.5); }

    ///
    /// Set everyting to default values, here it means that everything
    /// is set to zero except the boundary values that are kept.
    ///
    void reset(void);

    ///
    /// @return The width of the bins in the histogram.
    ///
    inline double spacing(void) const { return (xmax_-xmin_)/nof_bins(); }

    ///
    /// @return The histogram upper boundary.
    ///
    /// @note The upper boundary value is outside the histogram.
    ///
    inline double xmax(void) const { return xmax_; }

    ///
    /// @return The histogram lower boundary.
    ///
    /// @note The lower boundary value is inside the histogram.
    ///
    inline double xmin(void) const { return xmin_; }

    ///
    /// @return The count of bin \a k in the histogram.
    ///
    inline double operator[](size_t k) const { return histogram_[k]; } 

    ///
    /// The assignment operator
    ///
    const Histogram& operator=(const Histogram&);

  private:
    // Returns zero if outside boundaries
    inline size_t bin(double d)
    { return (((d<xmin_) || (d>xmax_)) ? 0 :
              static_cast<size_t>(floor((d-xmin_)/spacing() ))); }

    std::vector<double> histogram_;
    double xmax_;
    double xmin_;
    statistics::AveragerWeighted sum_all_;      // average of all data
    statistics::AveragerWeighted sum_histogram_;// average of data in histogram
  };

///
/// The Histogram output operator
///
std::ostream& operator<<(std::ostream& s,const Histogram&);

}} // of namespace statistics and namespace theplu

#endif