00001 #ifndef _MISC_MATH_H
00002 #define _MISC_MATH_H
00003
00004 #include <numeric>
00005 #include <cmath>
00006 #include <algorithm>
00007 #include <functional>
00008
00009 namespace osl
00010 {
00011 namespace misc
00012 {
00016 template<unsigned int N, class T>
00017 T nthPower(T x)
00018 {
00019 if (N==0)
00020 return 1;
00021
00022 const T temp = nthPower<N/2>(x);
00023 if (N%2 == 0)
00024 return temp * temp;
00025 else
00026 return temp * temp * x;
00027 }
00028
00029 template<class T, int N>
00030 struct SumDiffNthPower
00031 {
00032 T mean;
00033 SumDiffNthPower(T mean) : mean(mean) {}
00034 T operator()(T sum, T current)
00035 {
00036 return sum + nthPower<N>(current - mean);
00037 }
00038 };
00039
00040 template<class T, int N, class Iter_T>
00041 T nthMoment(Iter_T first, Iter_T last, T mean)
00042 {
00043 const size_t cnt = distance(first, last);
00044 return accumulate(first, last, T(), SumDiffNthPower<T, N>(mean))/cnt;
00045 }
00046
00047 template<class T, class Iter_T>
00048 T computeVariance(Iter_T first, Iter_T last, T mean)
00049 {
00050 return nthMoment<T, 2>(first, last, mean);
00051 }
00052
00053 template<class T, class Iter_T>
00054 T computeStdDev(Iter_T first, Iter_T last, T mean)
00055 {
00056 return sqrt(computeVariance(first, last, mean));
00057 }
00058
00059 template<class T, class Iter_T>
00060 T computeSkew(Iter_T first, Iter_T last, T mean)
00061 {
00062 const T m3 = nthMoment<T, 3>(first, last, mean);
00063 const T m2 = nthMoment<T, 2>(first, last, mean);
00064 return m3 / (m2 * sqrt(m2));
00065 }
00066
00067 template<class T, class Iter_T>
00068 T computeKurtosisExcess(Iter_T first, Iter_T last, T mean)
00069 {
00070 const T m4 = nthMoment<T, 4>(first, last, mean);
00071 const T m2 = nthMoment<T, 2>(first, last, mean);
00072 return m4 / (m2 * m2) - 3;
00073 }
00074
00075 template<class T, class Iter_T>
00076 void computeStats(Iter_T first, Iter_T last,
00077 T& sum, T& mean, T&var, T&std_dev, T& skew, T& kurt)
00078 {
00079 const size_t cnt = distance(first, last);
00080 sum = accumulate(first, last, T());
00081 mean = sum / cnt;
00082 var = computeVariance(first, last, mean);
00083 std_dev = sqrt(var);
00084 skew = computeSkew(first, last, mean);
00085 kurt = computeKurtosisExcess(first, last, mean);
00086 }
00087 }
00088 }
00089 #endif
00090
00091
00092
00093