From 512a34161c2d9f2560fe38821d63059c805cda8e Mon Sep 17 00:00:00 2001 From: Vanya Belyaev Date: Sun, 7 Apr 2024 16:01:56 +0200 Subject: [PATCH] 1. Add functions for the 1st, 2nd3rd and 4th cumulants --- ReleaseNotes/release_notes.md | 1 + ostap/stats/moment.py | 137 ++++++++++++++++++++++++++++++++- source/include/Ostap/Combine.h | 48 ++++++++++-- source/include/Ostap/Moments.h | 130 ++++++++++++++++++++++++++----- source/src/Combine.cpp | 53 ++++++++++++- 5 files changed, 340 insertions(+), 29 deletions(-) diff --git a/ReleaseNotes/release_notes.md b/ReleaseNotes/release_notes.md index 3e9bcf9a..205ad8e2 100644 --- a/ReleaseNotes/release_notes.md +++ b/ReleaseNotes/release_notes.md @@ -1,4 +1,5 @@ ## New features + 1. Add functions for the 1st, 2nd3rd and 4th cumulants ## Backward incompatible diff --git a/ostap/stats/moment.py b/ostap/stats/moment.py index 668a7c45..66339301 100644 --- a/ostap/stats/moment.py +++ b/ostap/stats/moment.py @@ -107,7 +107,7 @@ def _om_skewness ( obj ) : # m = ... # v = m.kurtosis () # @endcode -def _om_kurtosis( obj ) : +def _om_kurtosis ( obj ) : """Get an excess kurtosis for the moment-counter >>> m = ... >>> v = m.kurtosis () @@ -119,6 +119,93 @@ def _om_kurtosis( obj ) : return m4 / ( m2 * m2 ) - 3.0 return Ostap.Math.Moments.kurtosis ( obj ) +# ============================================================================= +## Get 1st cumulant, well i tis equasl to the mean +# @code +# m = ... +# v = m.cumulant1() +# @endcode +def _om_cumulant_1st( obj ) : + """Get 1st cumulant, well it is equal to the mean + + >>> m = ... + >>> v = m.cumulant_1st() + """ + return obj.mean() + +# ============================================================================= +## Get 2nd cumulant, well it is equal to the variance +# @code +# m = ... +# v = m.cumulant2() +# @endcode +def _om_cumulant_2nd ( obj ) : + """Get 3nd cumulant, well it is equal to the variance + + >>> m = ... + >>> v = m.cumulant2() + """ + assert 2 <= obj.order , 'cumulant_2nd: the order must be >=2!' + ## + if root_info < ( 6 , 18 ) : + ## get the unbiased 2nd moment + return obj.variance () + ## + return Ostap.Math.Moments.cumulant_2nd ( obj ) + +# ============================================================================= +## Get 3nd cumulant, well it is equal to the 3rd central moment +# @code +# m = ... +# v = m.cumulant3() +# @endcode +def _om_cumulant_3rd ( obj ) : + """Get 3rd cumulant, well it is equal to the 3rd central moment + + >>> m = ... + >>> v = m.cumulant_3rd() + """ + assert 3 <= obj.order , 'cumulant_3rd: the order must be >=3!' + ## + if root_info < ( 6 , 18 ) : + ## get the unbiased 3rd moment + unb = m.unbiased_3rd () + if obj.order() < 6 : return unb + ## + ## get moment + val = m.central_moment ( 3 ) + ## get the correctedalue + if isinstance ( val , VE ) : return VE ( unb , val.cov2() ) + ## + return Ostap.Math.Moments.cumulant_3rd ( obj ) + + +# ============================================================================= +## Get 3th cumulant, +# @code +# m = ... +# v = m.cumulant3() +# @endcode +def _om_cumulant_4th ( obj ) : + """Get 34th cumulant, well it is equal to the 3rd central moment + + >>> m = ... + >>> v = m.cumulant_4th() + """ + assert 4 <= obj.order , 'cumulant_4th: the order must be >=3!' + ## + if root_info < ( 6 , 18 ) : + ## get the unbiased 3rd moment + m4 = obj.moment ( 4 ) + m2 = obj.moment ( 2 ) + n = obj.size() + vv = m4 * ( n - 1 ) - 3 * m2 * m2 * ( n - 1 ) + ## + return vv * n * n * 1.0 / ( ( n - 1 ) * ( n -2 ) * ( n -3 ) ) + ## + return Ostap.Math.Moments.cumulant_4th ( obj ) + + # ============================================================================= ## get unbiased 2nd order moment from the moment-counter # @code @@ -184,6 +271,7 @@ def _om_u5th ( obj ) : return Ostap.Math.Moments.unbiased_5th ( obj ) + pos_infinity = float('+Inf') neg_infinity = float('-Inf') @@ -410,6 +498,47 @@ def _om_table ( obj , title = '' , prefix = '' , standard = False ) : row = "M[5](unb)" , '' if not n else '[10^%+d]' % n , field rows.append ( row ) + if 1 <= order and 1 <= size and obj.ok () and hasattr ( obj , 'cumulant_1st' ) : + + v = obj.cumulant_1st () + vv = float ( v ) + if isfinite ( vv ) and IM != float ( v ) : + if isinstance ( v , VE ) : field , n = pretty_ve ( v ) + else : field , n = pretty_float ( v ) + row = "K[1](unb)" , '' if not n else '[10^%+d]' % n , field + rows.append ( row ) + + if 2 <= order and 2 <= size and obj.ok () and hasattr ( obj , 'cumulant_2nd' ) : + + v = obj.cumulant_2nd () + vv = float ( v ) + if isfinite ( vv ) and IM != float ( v ) : + if isinstance ( v , VE ) : field , n = pretty_ve ( v ) + else : field , n = pretty_float ( v ) + row = "K[2](unb)" , '' if not n else '[10^%+d]' % n , field + rows.append ( row ) + + if 3 <= order and 3 <= size and obj.ok () and hasattr ( obj , 'cumulant_3rd' ) : + + v = obj.cumulant_2nd () + vv = float ( v ) + if isfinite ( vv ) and IM != float ( v ) : + if isinstance ( v , VE ) : field , n = pretty_ve ( v ) + else : field , n = pretty_float ( v ) + row = "K[3](unb)" , '' if not n else '[10^%+d]' % n , field + rows.append ( row ) + + if 4 <= order and 4 <= size and obj.ok () and hasattr ( obj , 'cumulant_4th' ) : + + v = obj.cumulant_4th () + vv = float ( v ) + if isfinite ( vv ) and IM != float ( v ) : + if isinstance ( v , VE ) : field , n = pretty_ve ( v ) + else : field , n = pretty_float ( v ) + row = "K[4](unb)" , '' if not n else '[10^%+d]' % n , field + rows.append ( row ) + + fmt = '[%d]' if 10 <= order < 100 : fmt = '[%02d]' elif 100 <= order < 1000 : fmt = '[%03d]' @@ -462,6 +591,12 @@ def _om_table ( obj , title = '' , prefix = '' , standard = False ) : Ostap.Math.Moment.central_moment = _om_cm2 Ostap.Math.Moment.table = _om_table +Ostap.Math.Moment.cumulant_1st = _om_cumulant_1st ## 1st cumulant is a mens +Ostap.Math.Moment.cumulant_2nd = _om_cumulant_2nd ## 2nd cumulant is a variance +Ostap.Math.Moment.cumulant_3rd = _om_cumulant_3rd ## 3nd cumulant is a 3rd central moment +Ostap.Math.Moment.cumulant_4th = _om_cumulant_4th ## 4th cumulant is different + + Ostap.Math.WMoment.mean = _om_mean Ostap.Math.WMoment.rms = _om_rms Ostap.Math.WMoment.variance = _om_variance diff --git a/source/include/Ostap/Combine.h b/source/include/Ostap/Combine.h index 504e07cb..29751e49 100644 --- a/source/include/Ostap/Combine.h +++ b/source/include/Ostap/Combine.h @@ -84,7 +84,7 @@ namespace Ostap { if ( Ostap::Math::inverse ( m_cov2 , m_vxi ) ) { - Ostap::throwException ( "Covariance matrix is not innvertile!" , + Ostap::throwException ( "Covariance matrix is not innvertible!" , "Ostap::Math::Combine<>" , 730 ) ; } const Data& vone = this->units() ; @@ -131,7 +131,7 @@ namespace Ostap // ====================================================================== public: // ====================================================================== - /// the main method: get a combined value using the calculated weights + /// the main method: get a combined value using the calculated weights Ostap::Math::ValueWithError result () const { const double r = ROOT::Math::Dot ( m_data , m_w ) ; @@ -198,9 +198,47 @@ namespace Ostap * @date 2015-09-28 */ Ostap::Math::ValueWithError - combine ( const double x , - const double y , - const Ostap::SymMatrix2x2& cov ) ; + combine + ( const double x , + const double y , + const Ostap::SymMatrix2x2& cov ) ; + // ======================================================================== + /** combine three measurements x, y + * and z with covarinace matrix cov + * @param x (INPUT) the first measurement + * @param y (INPUT) the second measurement + * @param z (INPUT) the third measurement + * @param cov (INPUT) covariance matrix + * @return combined result + * @author Vanya BELYAEV Ivan.Belyaev@itep.ru + * @date 2015-09-28 + */ + Ostap::Math::ValueWithError + combine + ( const double x , + const double y , + const double z , + const Ostap::SymMatrix3x3& cov ) ; + // ======================================================================== + /** combine four measurements x, y, + * z and w + * with covarinace matrix cov + * @param x (INPUT) the first measurement + * @param y (INPUT) the second measurement + * @param z (INPUT) the third measurement + * @param w (INPUT) the fourth measurement + * @param cov (INPUT) covariance matrix + * @return combined result + * @author Vanya BELYAEV Ivan.Belyaev@itep.ru + * @date 2015-09-28 + */ + Ostap::Math::ValueWithError + combine + ( const double x , + const double y , + const double z , + const double w , + const Ostap::SymMatrix4x4& cov ) ; // ======================================================================== /** combine two measurements x1 and x2 * using correlation coefficient rho: \f$-1\le\rho\le1\f$ diff --git a/source/include/Ostap/Moments.h b/source/include/Ostap/Moments.h index 2978432e..c91135af 100644 --- a/source/include/Ostap/Moments.h +++ b/source/include/Ostap/Moments.h @@ -1235,7 +1235,7 @@ namespace Ostap static inline double unbiased_2nd ( const Moment_& m ) { const unsigned long long n = m.size() ; - return !m.ok() || n < 2 ? s_INVALID_MOMENT : m.template M_<2> () / ( n - 1 ) ; + return !m.ok() || n < 2 ? s_INVALID_MOMENT : m.template M_<2> () / ( n - 1 ) ; } // ====================================================================== /** get the unbiased estimator for the 3rd order moment: @@ -1248,7 +1248,7 @@ namespace Ostap static inline double unbiased_3rd ( const Moment_& m ) { const unsigned long long n = m.size() ; - return !m.ok() || n < 3 ? s_INVALID_MOMENT : m.template M_<3> * n / ( ( n - 1.0L ) * ( n - 2.0L ) ) ; + return !m.ok() || n < 3 ? s_INVALID_MOMENT : m.template M_<3> * n / ( ( n - 1.0L ) * ( n - 2.0L ) ) ; } // ====================================================================== /** get the unbiased estimator for the 4th order moment: @@ -1266,7 +1266,8 @@ namespace Ostap static inline double unbiased_4th ( const Moment_& m ) { const unsigned long long n = m.size() ; - if ( !m.ok() || 4 > n ) { return s_INVALID_MOMENT ; } + if ( !m.ok() || (n < 4 ) ) { return s_INVALID_MOMENT ; } + // const long double m4 = m.template M_<4> / n ; const long double m2 = m.template M_<2> / n ; // @@ -1285,7 +1286,8 @@ namespace Ostap static inline double unbiased_5th ( const Moment_& m ) { const unsigned long long n = m.size() ; - if ( !m.ok() || 5 > n ) { return s_INVALID_MOMENT ; } + if ( !m.ok() || ( n < 5 ) ) { return s_INVALID_MOMENT ; } + // const long double m5 = m.template M_<5> () / n ; const long double m2 = m.template M_<2> () / n ; const long double m3 = m.template M_<3> () / n ; @@ -1330,7 +1332,7 @@ namespace Ostap // if ( n < 4 ) { return m2 ; } // ATTENTION! // - const double m4 = m.template M_<4> () / n ; + const double m4 = m.template M_<4> () / n ; // const double cov2 = ( m4 - m2 * m2 * ( n - 3.0L ) / ( n - 1.0L ) ) / n ; // @@ -1346,7 +1348,8 @@ namespace Ostap // const double m3 = unbiased_3rd ( m ) ; const double m2 = m.template M_<2> () / n ; - const double skew = m3 / std::pow ( m2 , 3.0/2 ) ; + const double skew = m3 / std::pow ( m2 , 3.0/2 ) ; + // const double cov2 = 6.0L * n * ( n - 1 ) / ( ( n - 2.0L ) * ( n + 1.0L ) * ( n + 3.0L ) ) ; return VE ( skew , cov2 ) ; } @@ -1357,16 +1360,17 @@ namespace Ostap { const unsigned long long n = m.size() ; if ( !m.ok() || n < 4 ) { return VE ( s_INVALID_MOMENT , -1 ) ;} - const double m4 = unbiased_4th ( m ) ; - const double m2 = m.template M_<2> () / n ; + const double m4 = unbiased_4th ( m ) ; + const double m2 = m.template M_<2> () / n ; // - const double k = m4 / ( m2 * m2 ) - 3 ; - double cov2 = 6.0L * n * ( n - 1 ) / ( ( n - 2.0L ) * ( n + 1.0L ) * ( n + 3.0L ) ) ; + const double k = m4 / ( m2 * m2 ) - 3 ; + double cov2 = 6.0L * n * ( n - 1 ) / ( ( n - 2.0L ) * ( n + 1.0L ) * ( n + 3.0L ) ) ; cov2 *= 4.0L * ( n * 1.0L * n -1 ) / ( ( n - 3.0L ) * ( n + 5.0L ) ) ; // return VE ( k , cov2 ) ; - } - + } + // ====================================================================== + // ====================================================================== /// get the central moment of order \f$ N \f$ template =2*K),int>::type = 0 > static inline VE central_moment ( const Moment_& m ) { return ( !m.ok() || m.size() < K ) ? VE(s_INVALID_MOMENT,-1) : m.template moment_ () ; } + // ====================================================================== // ====================================================================== - /// get the standartized moment of order 1 + /// get the standartized central moment of order 1 template ::type = 1 > static inline double std_moment ( const Moment_& /* m */ ) { return 1 ; } - /// get the central moment of order \f$ N \f$ + /// get the standartised central moment of order \f$ K \f$ template ::type = 1 > static inline double std_moment ( const Moment_& /* m */ ) { return 0 ; } - /// get the central moment of order \f$ N \f$ + /// get the standartised central moment of order \f$ N \f$ template =2)&&(K<=N),int>::type = 1 > static inline double std_moment ( const Moment_& m ) { return ( !m.ok() || m.size() < K ) ? s_INVALID_MOMENT : m.template moment_ () ; } - + // ====================================================================== + // ====================================================================== - // Weighted + /// get the first cumulant, well, it is actually the first central moment + template ::type = 1 > + static inline double cumulant_1st ( const Moment_& m ) + { return ( !m.ok() || ( m.size() < 1 ) ) ? s_INVALID_MOMENT : m.mu () ; } + /// get the first cumulant, well, it is actually the first central moment + template 1),int>::type = 1 > + static inline VE cumulant_1st ( const Moment_& m ) + { + if ( !m.ok() || m.size() < 1 ) { return s_INVALID_MOMENT ; } + // + const double m1 = m.mu() ; + if ( m.size () < 2 ) { return m1 ; } + // + const double m2 = m.template moment_<2>() ; + const double c2 = m2 / m.size() ; + // + return VE ( m1 , c2 ) ; + } + /// get the second unbiased cumulant, well it is actually the second unbiased central moment + template =2)&&(N<=3),int>::type = 1 > + static inline double cumulant_2nd ( const Moment_& m ) + { return ( !m.ok() || ( m.size() < 2 ) ) ? s_INVALID_MOMENT : m.size() * m.template M_<2> () / ( m.size() - 1 ) ; } + /// get the second unbiased cumulant, well it is actually the second unbiased central moment + template =4),int>::type = 1 > + static inline VE cumulant_2nd ( const Moment_& m ) + { + if ( !m.ok() || ( m.size() < 2 ) ) { return s_INVALID_MOMENT ;} + // + const auto n = m.size() ; + // + const double k2 = m.template M_<2> () / ( n - 1 ) ; // unbinased estiimate + if ( m.size() < 4 ) { return k2 ; } + // + const double m2 = m.template M_<2> () / n ; + const double m4 = m.template M_<4> () / n ; + // + const double k4 = + ( ( n + 1 ) * m4 - 3 * m2 * m2 * ( n - 1 ) ) * n * n + / ( ( n -1 ) * ( n - 2 ) * ( n - 3 ) ) ; + // unbinased variance + const double c2 = ( 2 * k2 * k2 * n + ( n - 1 ) * k4 ) / ( n * ( n + 1 ) ) ; + // + return VE ( k2 , c2 ) ; + } + /// get the third unbiased cumulant, well it is actually the third unbiased central moment + template =3),int>::type = 1 > + static inline double cumulant_3rd ( const Moment_& m ) + { + if ( !m.ok() || ( m.size() < 3 ) ) { return s_INVALID_MOMENT ; } + // + const auto n = m.size() ; + const double m3 = m.template M_<3> () / n ; + // + return m3 * n * n / ( ( n - 1 ) * ( n - 2 ) ) ; + } + /// get the 4th unbiased cumulant + template =3),int>::type = 1 > + static inline double cumulant_4th ( const Moment_& m ) + { + if ( !m.ok() || ( m.size() < 4 ) ) { return s_INVALID_MOMENT ; } + // + const auto n = m.size() ; + // + const double m2 = m.template M_<2> () / n ; + const double m4 = m.template M_<4> () / n ; + // + const double k4 = ( ( n + 1 ) * m4 - 3 * m2 *m2 * ( n - 1 ) ) + / ( ( n - 1 ) * ( n - 2 ) * ( n -3 ) ) ; + // + return k4 ; + } // ====================================================================== + // ====================================================================== + // Weighted + // ====================================================================== + // ====================================================================== /// get the mean static inline double mean ( const WMoment_<1>& m ) { return m.mu () ; } @@ -1428,7 +1514,7 @@ namespace Ostap /// get the variance template N),int>::type = 0 > static inline double variance ( const WMoment_& m ) - { return !m.ok() || m.size() < 2 ? s_INVALID_MOMENT : m.template moment_<2>() ; } + { return !m.ok() || ( m.size() < 2 ) ? s_INVALID_MOMENT : m.template moment_<2>() ; } // ====================================================================== /// get the sample variance with uncertainty template 3),int>::type = 0 > @@ -1455,7 +1541,8 @@ namespace Ostap template =3),int>::type = 0 > static inline VE skewness ( const WMoment_& m ) { - if ( !m.ok() || m.size() < 3 ) { return VE ( s_INVALID_MOMENT , -1 ) ; } + if ( !m.ok() || ( m.size() < 3 ) ) { return VE ( s_INVALID_MOMENT , -1 ) ; } + // const auto n = m.nEff () ; const double m3 = m.template M_<3> () / m.w () ; const double m2 = m.template M_<2> () / m.w () ; @@ -1469,7 +1556,8 @@ namespace Ostap template 3),int>::type = 0 > static inline VE kurtosis ( const WMoment_& m ) { - if ( !m.ok() || m.size() < 4 ) { return VE ( s_INVALID_MOMENT , -1 ) ; } + if ( !m.ok() || ( m.size() < 4 ) ) { return VE ( s_INVALID_MOMENT , -1 ) ; } + // const auto n = m.nEff () ; const double m4 = m.template M_<4> () / m.w() ; const double m2 = m.template M_<2> () / m.w() ; @@ -1479,8 +1567,8 @@ namespace Ostap // return 0 <= cov2 ? VE ( k , cov2 ) : VE ( k , 0.0 ) ; } + // ====================================================================== - // ====================================================================== /// get the central moment of order \f$ N \f$ template combiner ( data , cov ) ; + const Ostap::Math::Combine<2> combiner ( data , cov ) ; + return combiner.result() ; +} +// ============================================================================ +/* combine three measurements x, y + * and z with covarinace matrix cov + * @param x (INPUT) the first measurement + * @param y (INPUT) the second measurement + * @param z (INPUT) the third measurement + * @param cov (INPUT) covariance matrix + * @return combined result + * @author Vanya BELYAEV Ivan.Belyaev@itep.ru + * @date 2015-09-28 + */ +// ============================================================================ +Ostap::Math::ValueWithError +Ostap::Math::combine +( const double x , + const double y , + const double z , + const Ostap::SymMatrix3x3& cov ) +{ + Ostap::Vector3 data ( x , y , z ) ; + const Ostap::Math::Combine<3> combiner ( data , cov ) ; + return combiner.result() ; +} +// ============================================================================ +/* combine four measurements x, y, + * z and w + * with covarinace matrix cov + * @param x (INPUT) the first measurement + * @param y (INPUT) the second measurement + * @param z (INPUT) the third measurement + * @param w (INPUT) the fourth measurement + * @param cov (INPUT) covariance matrix + * @return combined result + * @author Vanya BELYAEV Ivan.Belyaev@itep.ru + * @date 2015-09-28 + */ +// ============================================================================ +Ostap::Math::ValueWithError +Ostap::Math::combine +( const double x , + const double y , + const double z , + const double w , + const Ostap::SymMatrix4x4& cov ) +{ + Ostap::Vector4 data ( x , y , z , w ) ; + const Ostap::Math::Combine<4> combiner ( data , cov ) ; return combiner.result() ; } // ============================================================================