---

more/phys/confidence_interval.h

Go to the documentation of this file.

//  Copyright (C) 2001--2002  Petter Urkedal (petter.urkedal@matfys.lth.se)
//
//  This file 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.
//
//  This file 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
//
//  As a special exception, you may use this file as part of a free
//  software library without restriction.  Specifically, if other files
//  instantiate templates or use macros or inline functions from this
//  file, or you compile this file and link it with other files to
//  produce an executable, this file does not by itself cause the
//  resulting executable to be covered by the GNU General Public
//  License.  This exception does not however invalidate any other
//  reasons why the executable file might be covered by the GNU General
//  Public License.
//
//  $Id: confidence_interval.h,v 1.1 2002/05/30 18:01:40 petter_urkedal Exp $


#ifndef MORE_PHYS_CONFIDENCE_INTERVAL_H
#define MORE_PHYS_CONFIDENCE_INTERVAL_H


#include <iosfwd>
#include <cmath>
#include <cstdlib>
#include <limits>
#include <more/math/math.h>
#include <more/io/syncstream.h>


namespace more {
namespace phys {

  /** Origin of a quantity.  Origins may be ored together, except for
   ** \c origin_theory which is 0 and \c origin_unknown which is the
   ** binary or of all the others.  */
  enum origin_t
  {
      /** Use for mathematical quantities or physical quantities which
       ** are exact by definition or known from a trusted theory (but
       ** not calculated from a model).  This origin is exclusive (can
       ** not be combined with the other origins). */
      origin_theory             = 0,

      /** The value is based on an experiment, and may usually
       ** be trusted within some small multiple of the given
       ** deviation. */
      origin_experiment         = 1,

      /** The value is based on a calculation and depends on
       ** the correctness of the model used. */
      origin_calculation        = 2,

      /** The value is based on systematics and may therefore
       ** have unexpected deviations. */
      origin_systematics        = 4,

      /** An origin different from the above.  This may have a
       ** specific meaning in client code or in parts of the code. */
      origin_other              = 8,

      /** Use if origin is unknown, uspecified, or other than
       ** the below.  This origin implies the other
       ** non-theoretic origins and is the least reliable. */
      origin_unknown            = 15
  };

  static unsigned int const origin_width = 4;

  inline origin_t operator|(origin_t x, origin_t y)
  { return origin_t(static_cast<int>(x) | static_cast<int>(y)); }
  inline origin_t operator&(origin_t x, origin_t y)
  { return origin_t(static_cast<int>(x) & static_cast<int>(y)); }
  inline origin_t operator^(origin_t x, origin_t y)
  { return origin_t(static_cast<int>(x) ^ static_cast<int>(y)); }
  inline origin_t operator~(origin_t x)
  { return origin_t(~ static_cast<int>(x)); }

  inline origin_t operator|=(origin_t& x, origin_t y) { return x = x | y; }
  inline origin_t operator&=(origin_t& x, origin_t y) { return x = x & y; }
  inline origin_t operator^=(origin_t& x, origin_t y) { return x = x ^ y; }


  /** \class confidence_interval confidence_interval.h \
   **     more/phys/confidence_interval.h
   **
   ** A type for storing a value and its deviation.
   ** \pre \c T is a floating point number. <tt>Moment > 0</tt>. */
  template<typename T>
    struct confidence_interval
    {
        typedef origin_t origin_type;
        static struct unknown_tag {} unknown;
        static struct infinite_tag {} infinite;

      private:
        enum pf_t
        {
            pf_unknown = 0,
            pf_finite = 1,
            pf_infinite = 2
        };
      public:
        /// Constuct an completely undefined value.
        confidence_interval()
            : m_origin(origin_unknown),
              m_pf_cent(pf_unknown),
              m_pf_s_mi(pf_unknown),
              m_pf_s_pl(pf_unknown),
              m_aux_flags(0),
              m_moment(0),
              m_cent(nan()),
              m_s_mi(nan()),
              m_s_pl(nan()) {}

        /// Construct an infinite value.
        confidence_interval(infinite_tag, origin_type origin = origin_unknown,
                            unsigned int moment = 2)
            : m_origin(origin),
              m_pf_cent(pf_infinite),
              m_pf_s_mi(pf_finite),
              m_pf_s_pl(pf_finite),
              m_aux_flags(0),
              m_moment(moment),
              m_cent(nan()),
              m_s_mi(0.0),
              m_s_pl(0.0) {}

        /** Construct an "approximately infinite" value.  Typically
         ** this may express a value which are unknown because it is
         ** too large too detect with present technology. */
        confidence_interval(infinite_tag, unknown_tag,
                            origin_type origin = origin_unknown,
                            unsigned int moment = 2)
            : m_origin(origin),
              m_pf_cent(pf_infinite),
              m_pf_s_mi(pf_unknown),
              m_pf_s_pl(pf_unknown),
              m_aux_flags(0),
              m_moment(moment),
              m_cent(nan()),
              m_s_mi(0.0),
              m_s_pl(0.0) {}

        /// Construct a symmetric confidence interval centered at \a x
        /// with standard deviation \a s.
        confidence_interval(T const& x, T const& s,
                            origin_type origin = origin_unknown,
                            unsigned int moment = 2)
            : m_origin(origin),
              m_pf_cent(pf_finite),
              m_pf_s_mi(pf_finite),
              m_pf_s_pl(pf_finite),
              m_aux_flags(0),
              m_moment(moment),
              m_cent(x),
              m_s_mi(s),
              m_s_pl(s) {}

        /// Construct an asymmetric confidence interval centered at \a
        /// x with deviations \a m_mi below and \a s_pl above.
        confidence_interval(T const& x, T const& s_mi, T const& s_pl,
                            origin_type origin = origin_unknown,
                            unsigned int moment = 2)
            : m_origin(origin),
              m_pf_cent(pf_finite),
              m_pf_s_mi(pf_finite),
              m_pf_s_pl(pf_finite),
              m_aux_flags(0),
              m_moment(moment),
              m_cent(x),
              m_s_mi(s_mi),
              m_s_pl(s_pl) {}

        /// Construct a lower bound \a x with deviation \a s_mi.
        confidence_interval(T const& x, T const& s_mi, infinite_tag,
                            origin_type origin = origin_unknown,
                            unsigned int moment = 2)
            : m_origin(origin),
              m_pf_cent(pf_finite),
              m_pf_s_mi(pf_finite),
              m_pf_s_pl(pf_infinite),
              m_aux_flags(0),
              m_moment(moment),
              m_cent(x),
              m_s_mi(s_mi),
              m_s_pl(inf()) {}

        /// Construct an upper bound \a x with deviation \a s_pl.
        confidence_interval(T const& x, infinite_tag, T const& s_pl,
                            origin_type origin = origin_unknown,
                            unsigned int moment = 2)
            : m_origin(origin),
              m_pf_cent(pf_finite),
              m_pf_s_mi(pf_infinite),
              m_pf_s_pl(pf_finite),
              m_aux_flags(0),
              m_moment(moment),
              m_cent(x),
              m_s_mi(nan()),
              m_s_pl(s_pl) {}

        /** Construct an approximate value \a x with unspecified
         ** deviation.  Approximate values are synced with 10 bits
         ** precition (œôøå 3 digits). */
        confidence_interval(T const& x, unknown_tag,
                            origin_type origin = origin_unknown)
            : m_origin(origin),
              m_pf_cent(pf_finite),
              m_pf_s_mi(pf_unknown),
              m_pf_s_pl(pf_unknown),
              m_aux_flags(0),
              m_moment(0),
              m_cent(x),
              m_s_mi(nan()),
              m_s_pl(nan()) {}

        /// Return generic information on the origin of the value.
        origin_type origin() const { return (origin_type)m_origin; }

        /// Set the generic origin of the value.
        void set_origin(origin_type o) { m_origin = o; }

        void be_uncertain() { m_aux_flags |= 1; }
        bool is_uncertain() const { return m_aux_flags & 1; }
        void be_signless() { m_aux_flags |= 2; }
        bool is_signless() const { return m_aux_flags & 2; }

//      void set_systematic() { m_origin |= origin_systematics; }
//      void set_calculated() { m_origin |= origin_calculation; }
//      bool is_systematic() const { return m_origin & origin_systematics; }
//      bool is_calculated() const { return m_origin & origin_calculation; }

        /// True if the value is defined is some way.
        bool is_known() const { return m_pf_cent != pf_unknown; }

        /// True if the valus is known and finite.
        bool is_finite() const { return m_pf_cent == pf_finite; }

        /// True if an infinite value.
        bool is_inf() const { return m_pf_cent == pf_infinite; }

        /// True if the lower deviation is given.
        bool dev_below_is_known() const { return m_pf_s_mi != pf_unknown; }

        /// True if the upper deviation is given.
        bool dev_above_is_known() const { return m_pf_s_pl != pf_unknown; }

        /// Deviation below is known and finite.
        bool dev_below_is_finite() const { return m_pf_s_mi == pf_finite; }

        /// Deviation above is known and finite.
        bool dev_above_is_finite() const { return m_pf_s_pl == pf_finite; }

        /// Deviation below is known and infinite.
        bool dev_below_is_inf() const { return m_pf_s_mi == pf_infinite; }

        /// Deviation above is known and infinite.
        bool dev_above_is_inf() const { return m_pf_s_pl == pf_infinite; }

        /// True if the deviation above and below are the same.
        bool is_symmetric() const
        {
            return m_pf_s_mi == m_pf_s_pl
                && (m_pf_s_mi != pf_finite || m_s_pl == m_s_mi);
        }

        /// True if exact within the numeric precision.
        bool is_exact() const
        {
            return dev_below_is_finite()
                && dev_above_is_finite()
                && m_s_mi == 0.0 && m_s_pl == 0.0;
        }

        /// Center of the confidence interval.
        T cent() const { return m_cent; }

        /// Minimum value in the confidence interval.
        T min() const { return m_cent - m_s_mi; }

        /// Maximum value in the confidence interval.
        T max() const { return m_cent + m_s_pl; }

        /// Return <tt>max(dev_above(), dev_below())</tt>.
        T dev() const
        {
            return m_s_mi > m_s_pl? m_s_mi : m_s_pl;
        }

        /// Return \c abs(cent())/dev().
        T rdev() const { return dev()/std::fabs(m_cent); }

        /// The deviation below the central value.
        T dev_below() const { return m_s_mi; }

        /// The deviation above the central value.
        T dev_above() const { return m_s_pl; }

        confidence_interval&
        operator*=(T x)
        {
            if (x < 0) {
                negate();
                x = -x;
            }
            m_cent *= x;
            m_s_pl *= x;
            m_s_mi *= x;
            return *this;
        }

        confidence_interval&
        operator/=(T x) { return (*this) *= T(1)/x; }

        confidence_interval&
        operator+=(confidence_interval const& rhs)
        {
            m_cent += rhs.m_cent;
            add_deviation(rhs);
            return *this;
        }

        confidence_interval&
        operator-=(confidence_interval const& rhs)
        {
            m_cent -= rhs.m_cent;
            add_deviation(rhs);
            return *this;
        }

        void negate()
        {
            m_cent = -m_cent;
            T tmp = m_s_pl;
            m_s_pl = m_s_mi;
            m_s_mi = tmp;
        }

        void sync(io::syncstream&);

    private:
        static T nan() { return std::numeric_limits<T>::quiet_NaN(); }
        static T inf() { return std::numeric_limits<T>::infinity(); }

        void add_deviation(confidence_interval const& rhs);

        unsigned int m_origin : origin_width;
        unsigned int m_pf_cent : 2;
        unsigned int m_pf_s_mi : 2;
        unsigned int m_pf_s_pl : 2;
        unsigned int m_aux_flags : 4;
        unsigned char m_moment;
        T m_cent;
        T m_s_mi;
        T m_s_pl;
    };


  template<typename T>
    inline confidence_interval<T>
    operator-(confidence_interval<T> x)
    {
        x.negate();
        return x;
    }
  template<typename T>
    inline confidence_interval<T>
    operator+(confidence_interval<T> const& x) { return x; }

  template<typename T>
    inline confidence_interval<T>
    operator*(confidence_interval<T> x, T const& y)
    {
        return x *= y;
    }
  template<typename T>
    inline confidence_interval<T>
    operator*(T const& x, confidence_interval<T> y)
    {
        return y *= x;
    }
  template<typename T>
    inline confidence_interval<T>
    operator/(confidence_interval<T> x, T const& y)
    {
        return x /= y;
    }
  template<typename T>
    inline confidence_interval<T>
    operator+(confidence_interval<T> x,
              confidence_interval<T> y)
    {
        return x += y;
    }
  template<typename T>
    inline confidence_interval<T>
    operator-(confidence_interval<T> x,
              confidence_interval<T> y)
    {
        return x -= y;
    }

  template<typename T>
    inline bool
    better_precision(confidence_interval<T> const& x,
                     confidence_interval<T> const& y)
    {
        if (x.is_indet_below() || x.is_indet_above())
            return y.is_indet_below() && y.is_indet_above();
        else if (y.is_indet_below() || y.is_indet_above())
            return true;
        else
            return x.dev() < y.dev();
    }

  template<typename Char, typename Traits, typename T>
    std::basic_ostream<Char, Traits>&
    operator<<(std::basic_ostream<Char, Traits>&,
               confidence_interval<T> const&);

  void set_origin_indicator(std::ostream&, origin_t ogn, char const* ind);

  void set_origin_indicator(std::ostream&, char const* pfx, char const* sfx,
                            char const* sep);

}}

#endif

---