Changeset 588

Timestamp:

Jun 21, 2006, 5:43:29 PM (17 years ago)

Author:

Markus Ringnér

Message:

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

Location:

trunk

Files:

• c++_tools/statistics/utility.h (modified) (4 diffs)
• test/statistics_test.cc (modified) (2 diffs)

trunk/c++_tools/statistics/utility.h

@@ -11 +11 @@
 #include <vector>
 
+#include <gsl/gsl_statistics_double.h>
+
 namespace theplu {
 namespace statistics {  
@@ -21 +23 @@
 
   ///
-  /// Calculates the probabilty to get \a k or smaller from a
+  /// Calculates the probability to get \a k or smaller from a
   /// hypergeometric distribution with parameters \a n1 \a n2 \a
   /// t. Hypergeomtric situation you get in the following situation:
@@ -31 +33 @@
   ///
   double cdf_hypergeometric_P(u_int k, u_int n1, u_int n2, u_int t);
+
+
+  ///
+  /// Computes the kurtosis of the data in a vector. The kurtosis
+  /// measures how sharply peaked a distribution is, relative to its
+  /// width. The kurtosis is normalized to zero for a gaussian
+  /// distribution.
+  ///
+  inline double kurtosis(const gslapi::vector& v)
+  {
+    const gsl_vector* gvp=v.gsl_vector_p();
+    return gsl_stats_kurtosis(gvp->data,gvp->stride,gvp->size);
+  }
 
 
@@ -110 +125 @@
                     const bool sorted=false);
 
+  ///
+  /// Computes the skewness of the data in a vector. The skewness
+  /// measures the asymmetry of the tails of a distribution.
+  ///
+  inline double skewness(const gslapi::vector& v)
+  {
+    const gsl_vector* gvp=v.gsl_vector_p();
+    return gsl_stats_skew(gvp->data,gvp->stride,gvp->size);
+  }
+  
+
 }} // of namespace statistics and namespace theplu
 

trunk/test/statistics_test.cc

@@ -7 +7 @@
 #include <cstdlib>
 #include <iostream>
-
+#include <cmath>
 
 int main()
@@ -18 +18 @@
     gsl_vec(i)=i;
   }
+
   double m=statistics::median(data);
   double m_gsl=statistics::median(gsl_vec);
   if (m!=4.5 || m!=m_gsl)
     return -1;
+
+  double tolerance=1e-10;
+  double skewness_gsl=statistics::skewness(gsl_vec);
+  if (fabs(skewness_gsl)>tolerance)
+    return -1;
+  double kurtosis_gsl=statistics::kurtosis(gsl_vec);
+  if (fabs(kurtosis_gsl+1.5616363636363637113)>tolerance) 
+    return -1;  
   return 0;
 }