BasicStatistics.h

Go to the documentation of this file.

// Copyright (c) 2002-2023, The OpenMS Team -- EKU Tuebingen, ETH Zurich, and FU Berlin
// SPDX-License-Identifier: BSD-3-Clause
//
// --------------------------------------------------------------------------
// $Maintainer: Timo Sachsenberg $
// $Authors: $
// --------------------------------------------------------------------------
 
#pragma once
 
#include <OpenMS/CONCEPT/Types.h>
#include <vector>
#include <ostream>
#include <cmath>
#include <numeric>
 
namespace OpenMS
{
  namespace Math
  {
 
    template <typename RealT = double>
    class
    BasicStatistics
    {
 
public:
 
      typedef RealT RealType;
 
      typedef std::vector<RealType> probability_container;
      typedef std::vector<RealType> coordinate_container;
 
      BasicStatistics() :
        mean_(0),
        variance_(0),
        sum_(0)
      {}
 
      BasicStatistics(BasicStatistics const & arg) :
        mean_(arg.mean_),
        variance_(arg.variance_),
        sum_(arg.sum_)
      {}
 
      BasicStatistics & operator=(BasicStatistics const & arg)
      {
        mean_      = arg.mean_;
        variance_  = arg.variance_;
        sum_       = arg.sum_;
        return *this;
      }
 
      void clear()
      {
        mean_ = 0;
        variance_ = 0;
        sum_ = 0;
      }
 
 
      template <typename ProbabilityIterator>
      void update(ProbabilityIterator probability_begin,
                  ProbabilityIterator const probability_end
                  )
      {
        clear();
        unsigned pos = 0;
        ProbabilityIterator iter = probability_begin;
 
        for (; iter != probability_end; ++iter, ++pos)
        {
          sum_  += *iter;
          mean_ += *iter * pos;
        }
        mean_ /= sum_;
 
        for (iter = probability_begin, pos = 0; iter != probability_end; ++iter, ++pos)
        {
          RealType diff = RealType(pos) - mean_;
          diff *= diff;
          variance_ += *iter * diff;
        }
        variance_ /= sum_;
 
        if (sum_ == 0 && (std::isnan(mean_) || std::isinf(mean_)) )
        {
          mean_ = 0;
          variance_ = 0;
        }
      }
 
 
      template <typename ProbabilityIterator, typename CoordinateIterator>
      void update(ProbabilityIterator const probability_begin,
                  ProbabilityIterator const probability_end,
                  CoordinateIterator  const coordinate_begin
                  )
      {
        clear();
        ProbabilityIterator prob_iter = probability_begin;
        CoordinateIterator  coord_iter = coordinate_begin;
 
        for (; prob_iter != probability_end; ++prob_iter, ++coord_iter)
        {
          sum_  += *prob_iter;
          mean_ += *prob_iter * *coord_iter;
        }
        mean_ /= sum_;
 
        for (prob_iter = probability_begin, coord_iter = coordinate_begin;
             prob_iter != probability_end;
             ++prob_iter, ++coord_iter
             )
        {
          RealType diff = *coord_iter - mean_;
          diff *= diff;
          variance_ += *prob_iter * diff;
        }
        variance_ /= sum_;
        return;
      }
 
      RealType mean() const { return mean_; }
      void setMean(RealType const & mean) { mean_ = mean; }
 
      RealType variance() const { return variance_; }
      void setVariance(RealType const & variance) { variance_ = variance; }
 
      RealType sum() const { return sum_; }
      void setSum(RealType const & sum) { sum_ = sum; }
 
 
      RealType normalDensity_sqrt2pi(RealType coordinate) const
      {
        coordinate -= mean();
        coordinate *= coordinate;
        return exp(-coordinate / RealType(2) / variance());
      }
 
      static RealType sqrt2pi() { return 2.50662827463100050240; }
 
      inline RealType normalDensity(RealType const coordinate) const
      {
        return normalDensity_sqrt2pi(coordinate) / sqrt2pi();
      }
 
      void normalApproximation(probability_container & probability)
      {
        normalApproximationHelper_(probability, probability.size());
        return;
      }
 
      void normalApproximation(probability_container & probability,
                               typename probability_container::size_type const size
                               )
      {
        probability.resize(size);
        normalApproximationHelper_(probability, size);
        return;
      }
 
      void normalApproximation(probability_container & probability,
                               coordinate_container  const & coordinate
                               )
      {
        probability.resize(coordinate.size());
        normalApproximationHelper_(probability, coordinate);
        return;
      }
 
      friend std::ostream & operator<<(std::ostream & os, BasicStatistics & arg)
      {
        os << "BasicStatistics:  mean=" << arg.mean() << "  variance=" << arg.variance() << "  sum=" << arg.sum();
        return os;
      }
 
protected:
 
 
 
      RealType mean_;
      RealType variance_;
      RealType sum_;
 
private:
 
 
 
      void normalApproximationHelper_(probability_container & probability,
                                      typename probability_container::size_type const size
                                      )
      {
        RealType gaussSum = 0;
        typename coordinate_container::size_type i;
 
        // precondition size == probability.size() is guaranteed by wrappers.
        for (i = 0; i < size; ++i)
        {
          gaussSum += normalDensity_sqrt2pi(RealType(i));
        }
 
        for (i = 0; i < size; ++i)
        {
          probability[i] = normalDensity_sqrt2pi(RealType(i)) / gaussSum * sum();
        }
        return;
      }
 
      void normalApproximationHelper_(probability_container & probability,
                                      coordinate_container const & coordinate
                                      )
      {
        RealType gaussSum = 0;
        typename coordinate_container::size_type i;
        typename coordinate_container::size_type const size = coordinate.size();
 
        for (i = 0; i < size; ++i)
        {
          gaussSum += normalDensity_sqrt2pi(coordinate[i]);
        }
 
        for (i = 0; i < size; ++i)
        {
          probability[i] = normalDensity_sqrt2pi(coordinate[i]) / gaussSum * sum();
        }
        return;
      }
 
 
    };
 
  }   // namespace Math
 
} // namespace OpenMS