Changeset 588


Ignore:
Timestamp:
Jun 21, 2006, 5:43:29 PM (15 years ago)
Author:
Markus Ringnér
Message:

Added statistics functions to calculate the higher moments: skewness and kurtosis using gsl functions

Location:
trunk
Files:
2 edited

Legend:

Unmodified
Added
Removed
  • trunk/c++_tools/statistics/utility.h

    r519 r588  
    1111#include <vector>
    1212
     13#include <gsl/gsl_statistics_double.h>
     14
    1315namespace theplu {
    1416namespace statistics { 
     
    2123
    2224  ///
    23   /// Calculates the probabilty to get \a k or smaller from a
     25  /// Calculates the probability to get \a k or smaller from a
    2426  /// hypergeometric distribution with parameters \a n1 \a n2 \a
    2527  /// t. Hypergeomtric situation you get in the following situation:
     
    3133  ///
    3234  double cdf_hypergeometric_P(u_int k, u_int n1, u_int n2, u_int t);
     35
     36
     37  ///
     38  /// Computes the kurtosis of the data in a vector. The kurtosis
     39  /// measures how sharply peaked a distribution is, relative to its
     40  /// width. The kurtosis is normalized to zero for a gaussian
     41  /// distribution.
     42  ///
     43  inline double kurtosis(const gslapi::vector& v)
     44  {
     45    const gsl_vector* gvp=v.gsl_vector_p();
     46    return gsl_stats_kurtosis(gvp->data,gvp->stride,gvp->size);
     47  }
    3348
    3449
     
    110125                    const bool sorted=false);
    111126
     127  ///
     128  /// Computes the skewness of the data in a vector. The skewness
     129  /// measures the asymmetry of the tails of a distribution.
     130  ///
     131  inline double skewness(const gslapi::vector& v)
     132  {
     133    const gsl_vector* gvp=v.gsl_vector_p();
     134    return gsl_stats_skew(gvp->data,gvp->stride,gvp->size);
     135  }
     136 
     137
    112138}} // of namespace statistics and namespace theplu
    113139
  • trunk/test/statistics_test.cc

    r511 r588  
    77#include <cstdlib>
    88#include <iostream>
    9 
     9#include <cmath>
    1010
    1111int main()
     
    1818    gsl_vec(i)=i;
    1919  }
     20
    2021  double m=statistics::median(data);
    2122  double m_gsl=statistics::median(gsl_vec);
    2223  if (m!=4.5 || m!=m_gsl)
    2324    return -1;
     25
     26  double tolerance=1e-10;
     27  double skewness_gsl=statistics::skewness(gsl_vec);
     28  if (fabs(skewness_gsl)>tolerance)
     29    return -1;
     30  double kurtosis_gsl=statistics::kurtosis(gsl_vec);
     31  if (fabs(kurtosis_gsl+1.5616363636363637113)>tolerance)
     32    return -1; 
    2433  return 0;
    2534}
Note: See TracChangeset for help on using the changeset viewer.