11.5. Вычисление дисперсии, стандартного отклонения и других статистических функций

11.5. Вычисление дисперсии, стандартного отклонения и других статистических функций

Проблема

Требуется рассчитать значение одной или нескольких обычных статистических функций, например дисперсии (variance), стандартного отклонения (standard deviation), коэффициента асимметрии (skew) и эксцесса (kurtosis) для последовательности чисел.

Решение

Функцию accumulate из заголовочного файла <numeric> можно использовать для расчета многих статистических параметров, а не только для суммирования пользовательских объектов-функций. Пример 11.9 показывает, как можно вычислить значения некоторых важных статистические функций при помощи accumulate.

Пример 11.9. Статистические функции

#include <numeric>

#include <cmath>

#include <algorithm>

#include <functional>

#include <vector>

#include <iostream>

using namespace std;

template<int N, class T>

T nthPnwer(T x) {

 T ret = x;

 for (int i=1; i < N; ++i) {

  ret *= x;

 }

 return ret;

}

template<class T, int N>

struct SumDiffNthPower {

 SumDiffNthPower(T x) : mean_(x) {};

 T operator()(T sum, T current) {

  return sum + nthPower<N>(current - mean_);

 }

 T mean_;

};

template<class T, int N, class Iter_T>

T nthMoment(Iter_T first, Iter_T last, T mean) {

 size_t cnt = distance(first, last);

 return accumulate(first, last, T(), SumDiffNthPower<T, N>(mean)) / cnt;

}

template<class T, class Iter_T>

T computeVariance(Iter_T first, Iter_T last, T mean) {

 return nthMoment<T, 2>(first, last, mean);

}

template<class T, class Iter_T>

T computeStdDev(Iter_T first, Iter_T last, T mean) {

 return sqrt(computeVariance(first, last, mean));

}

template<class T, class Iter_T>

T computeSkew(Iter_T begin, Iter_T end, T mean) {

 T m3 = nthMoment<T, 3>(begin, end, mean);

 T m2 = nthMoment<T, 2>(begin, end, mean);

 return m3 / (m2 * sqrt(m2));

}

template<class T, class Iter_T>

T computeKurtosisExcess(Iter_T begin, Iter_T end, T mean) {

 T m4 = nthMoment<T, 4>(begin, end, mean);

 T m2 = nthMoment<T, 2>(begin, end, mean);

 return m4 / (m2 * m2) - 3;

}

template<class T, class Iter_T>

void computeStats(Iter_T first, Iter_T last, T& sum, T& mean,

 T& var, T& std_dev, T& skew, T& kurt) {

 size_t cnt = distance(first, last);

 sum = accumulate(first, last, T());

 mean = sum / cnt;

 var = computeVariance(first, last, mean);

 std_dev = sort(var);

 skew = computeSkew(first, last, mean);

 kurt = computeKurtosisExcess(first, last, mean);

}

int main() {

 vector<int> v;

 v.push_back(2);

 v.push_back(4);

 v.push_back(8);

 v.push_back(10);

 v.push_back(99);

 v.push_back(1);

 double sum, mean, var, dev, skew, kurt;

 computeStats(v.begin(), v.end(), sum, mean, var, dev, skew, kurt);

 cout << "count = " << v.size() << " ";

 cout << "sum = " << sum << " ";

 cout << "mean = " << mean << " ";

 cout << "variance = " << var << " ";

 cout << "standard deviation = " << dev << " ";

 cout << "skew = " << skew << " ";

 cout << "kurtosis excess = " << kurt << " ";

 cout << endl;

}

Программа примера 11.9 выдает следующий результат

count = 6

sum = 124

mean = 20.6667

variance = 1237.22

standard deviation = 35.1742

skew = 1.75664

kurtosis excess = 1.14171

Обсуждение

Некоторые наиболее важные статистические функции (например, дисперсия, стандартное отклонение, коэффициент асимметрии и эксцесс) определяются исходя из нормализованных выборочных моментов. Статистические функции определяются немного по-разному в различных текстах. Здесь мы используем несмещенные определения статистических функций, которые сведены в табл. 11.1.

Табл. 11.1. Определения статистических функций

Проще всего программировать статистические функции, определяя их через моменты. Поскольку используется несколько различных моментов, каждый из которых характеризуется целочисленной константой, я передаю эту константу как параметр шаблона. Это в целом позволяет компилятору генерировать более эффективный программный код, потому что это целочисленное значение известно на этапе компиляции.

Функция момента определяется при помощи математического оператора суммы. Во всех случаях, когда речь идет об этом операторе, следует иметь в виду функцию accumulate, определенную в заголовочном файле <numeric>. Существует две разновидности функции accumulate: одна подсчитывает сумму, используя operator+, а другая использует функтор суммирования, который вы должны предоставить. Ваш функтор суммирования будет принимать значение накопленной суммы и значение конкретного элемента последовательности.

Пример 11.10 иллюстрирует работу функции accumulate, показывая, как предоставленный пользователем функтор вызывается для каждого элемента последовательности.

Пример 11.10. Пример реализации функции accumulate

template<class Iter_T, class Value_T, class BinOp_T>

Iter_T accumulate(Iter_T begin, Iter_T end, Value_T value, BinOp_T op) {

 while (begin != end) {

  value = op(value, *begin++)

 }

 return value;

}