オンライン版各種平均テンプレートライブラリ(後述するオンライン版各種統計量テンプレートライブラリ statistics_online.hpp などを含みます)の averages_online.hpp を作成しました。これは平均値を一つ前の平均値から求めるもので、リアルタイムで効率的に各種平均を求めるのに適しています。まず、通常どおりに各種平均を求めるための式は以下の通りです。
相加平均(arithmetic average):
\[\mu_A = \dfrac1n\sum_{i=1}^n x_i\]
相乗平均(geometric average):
\[\mu_G = \sqrt[n]{\prod_{i=1}^n x_i} = \exp\left( \dfrac1n\sum_{i=1}^n \log x_i \right)\]
調和平均(harmonic average):
\[\mu_H = \dfrac n{\sum_{i=1}^n \dfrac1{x_i}}\]
一般化平均(generalized average):
\[\mu_{G(m)} = \sqrt[m]{\dfrac1n\sum_{i=1}^n x_i^m}\]
任意の一般化平均(arbitrary generalized average):
\[\mu_{f} = f^{-1}\left(\dfrac1n\sum_{i=1}^n f(x_i)\right)\]
以上より、一つ前の平均値から求めるには、以下の差分方程式に依ります。
相加平均の差分方程式(arithmetic average online):
\[\mu_A{(n)} = \dfrac1n\left\{ (n-1)\mu_A{(n-1)} + x_n \right\}\]
相乗平均の差分方程式(geometric average online):
\[\mu_G{(n)} = \sqrt[n]{\mu_G{(n-1)}^{n-1} x_n} = \exp\left( \dfrac1n\left\{ (n-1)\log(\mu_G{(n-1)}) + \log x_n \right\} \right)\]
調和平均の差分方程式(harmonic average online):
\[\mu_H{(n)} = \dfrac n{\dfrac{n-1}{\mu_H{(n-1)}} + \dfrac1{x_n}}\]
一般化平均の差分方程式(generalized average online):
\[\mu_{G(m)}{(n)} = \sqrt[m]{\dfrac1n\left\{ (n-1)\mu_{G(m)}{(n-1)}^m + x_n^m \right\}}\]
任意の一般化平均の差分方程式(arbitrary generalized average online):
\[\mu_{f}{(n)} = f^{-1}\left({\dfrac1n\left\{ (n-1)f\left(\mu_{f}{(n-1)}\right) + f(x_n) \right\}}\right)\]
また、もともと差分方程式で実装されるべき以下の平均もあります。
指数平滑平均(exponential smoothed average)の差分方程式:
\[\mu_E{(n)} = (1 - \alpha)\mu_E{(n-1)} + \alpha x_n\]
さらに、各種の移動平均は、幅\(T\)の記憶域から除去される最古の\(x_{n-T}\)と追加される最新の\(x_n\)および一つ前の移動平均値から求めることが出来、以下の差分方程式に依ります。
相加移動平均の差分方程式(moving arithmetic online):
\[\mu_{MA}{(n)} = \mu_{MA}{(n-1)} - \dfrac1Tx_{n-T} + \dfrac1Tx_n\]
相乗移動平均の差分方程式(moving geometric online):
\[\mu_{MG}{(n)} = \sqrt[T]{\mu_{MG}{(n-1)}^T\dfrac{x_n}{x_{n-T}}} = \exp\left(\log(\mu_{MG}{(n-1)}) - \dfrac1T\log(x_{n-T}) + \dfrac1T\log(x_n)\right)\]
調和移動平均の差分方程式(moving harmonic online):
\[\mu_{MH}{(n)} = \dfrac1{\dfrac1{\mu_{MH}{(n-1)}} - \dfrac1{Tx_{n-T}} + \dfrac1{Tx_n}}\]
一般化移動平均の差分方程式(moving generalized online):
\[\mu_{MG(m)}{(n)} = \sqrt[m]{\mu_{MG(m)}{(n-1)}^m - \dfrac1Tx_{n-T}^m + \dfrac1Tx_n^m}\]
任意の一般化移動平均(moving arbitrary generalized online)
\[\mu_{Mp}{(n)} = f^{-1}\left({f\left(\mu_{Mp}{(n-1)}\right) - \dfrac1Tf(x_{n-T}) + \dfrac1Tf(x_n)}\right)\]
ちなみに、その tarball には上記の検証の為に、任意のコンテナ(配列、vector, list 等々)で通常どおりに以下の平均を求める average.hpp も用意しました。
- 相加平均(arithmetic average)
- 相乗平均(geometric average)
- 調和平均(harmonic average)
- 一般化平均(generalized average)
- 任意の一般化平均(arbitrary generalized average)
- 指数平滑平均(exponential smoothed average)
これら非オンラインとオンラインとの自乗誤差は 100,000 データにおいて 10E-13 程度に収まっていることは確認できています。
使用例として、その tarball にあるコードを以下に示しておきます。
/*
averages.cc
Written by Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <vector>
#include <list>
#include <functional>
#include <cstdlib> // random(), etc
#include <cassert>
#include "averages.hpp"
#include "averages_online.hpp"
int main()
{
const double maximum_permissible_error(1e-14);
const char *result_names[] = {
"arithmetic mean",
"geometric mean",
"harmonic mean",
"generalized mean",
"exponential smoothed average",
};
std::vector<double> v(1000);
generic_average_online<double> *online_averages[] = {
new arithmetic_average_online<double>(),
new geometric_average_online<double>(),
new harmonic_average_online<double>(),
new generalized_average_online<double>(),
new exponential_smoothed_average_online<double>(),
};
srandomdev();
for (size_t i=0; i<v.size(); ++i) {
double x((double)random()/((1L<<31)-1));
v[i] = x;
for (size_t j=0; j<sizeof(online_averages)/sizeof(online_averages[0]); ++j)
online_averages[j]->average(x);
}
double results[] = {
arithmetic_average(v.begin(), v.end()),
geometric_average(v.begin(), v.end()),
harmonic_average(v.begin(), v.end()),
generalized_average(v.begin(), v.end()),
exponential_smoothed_average(v.begin(), v.end()),
};
for (size_t j=0; j<sizeof(online_averages)/sizeof(online_averages[0]); ++j) {
double error(sqrt(pow(online_averages[j]->get_value() - results[j], 2)));
std::cout << result_names[j] << ":\t" << online_averages[j]->get_value() << "\t" << results[j] << "\t" << error << std::endl;
assert(error < maximum_permissible_error);
}
const char *moving_result_names[] = {
"moving arithmetic average",
"moving geometric average",
"moving harmonic average",
"moving generalized average",
};
std::vector<double> m(v.size());
generic_average_online<double> *online_moving_averages[] = {
new moving_arithmetic_average_online<double>(v.size()),
new moving_geometric_average_online<double>(v.size()),
new moving_harmonic_average_online<double>(v.size()),
new moving_generalized_average_online<double>(v.size()),
};
for (size_t i=0; i<v.size()*10; ++i) {
double x((double)random()/((1L<<31)-1));
m[i % v.size()] = x;
for (size_t j=0; j<sizeof(online_moving_averages)/sizeof(online_moving_averages[0]); ++j)
online_moving_averages[j]->average(x);
if (i < v.size()) continue;
double moving_results[] = {
arithmetic_average(m.begin(), m.end()),
geometric_average(m.begin(), m.end()),
harmonic_average(m.begin(), m.end()),
generalized_average(m.begin(), m.end()),
};
for (size_t j=0; j<sizeof(online_moving_averages)/sizeof(online_moving_averages[0]); ++j) {
double error(sqrt(pow(online_moving_averages[j]->get_value() - moving_results[j], 2)));
if (i+1 == v.size()*10)
std::cout << moving_result_names[j] << ":\t" << online_moving_averages[j]->get_value() << "\t" << moving_results[j] << "\t" << error << std::endl;
assert(error < maximum_permissible_error);
}
}
return 0;
}
任意の一般化平均を追加したので、上例にて一致する箇所を入れ替える形で、その使用例を以下に示しておきます。
template <typename T>
struct identity_function : std::unary_function<T, T> {
T operator()(const T &x) const
{ return x; }
};
template <typename T>
struct reciprocal_function : std::unary_function<T, T> {
T operator()(const T &x) const
{ return T(1)/x; }
};
template <typename T>
struct square_function : std::unary_function<T, T> {
T operator()(const T &x) const
{ return x*x; }
};
#define USE_ARBITRARY_GENERALIZED_AVERAGE
#define USE_ARBITRARY_GENERALIZED_AVERAGE_ONLINE
generic_average_online<double> *online_averages[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE_ONLINE)
new arithmetic_average_online<double>(),
new geometric_average_online<double>(),
new harmonic_average_online<double>(),
new generalized_average_online<double>(),
#else
new arbitrary_generalized_average_online<double>(identity_function<double>(), identity_function<double>()),
new arbitrary_generalized_average_online<double>(log, exp),
new arbitrary_generalized_average_online<double>(reciprocal_function<double>(), reciprocal_function<double>()),
new arbitrary_generalized_average_online<double>(square_function<double>(), sqrt),
#endif
new exponential_smoothed_average_online<double>(),
};
double results[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE)
arithmetic_average(v.begin(), v.end()),
geometric_average(v.begin(), v.end()),
harmonic_average(v.begin(), v.end()),
generalized_average(v.begin(), v.end()),
#else
generalized_average(v.begin(), v.end(), identity_function<double>(), identity_function<double>()),
generalized_average(v.begin(), v.end(), log, exp),
generalized_average(v.begin(), v.end(), reciprocal_function<double>(), reciprocal_function<double>()),
generalized_average(v.begin(), v.end(), square_function<double>(), sqrt),
#endif
exponential_smoothed_average(v.begin(), v.end()),
};
generic_average_online<double> *online_moving_averages[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE_ONLINE)
new moving_arithmetic_average_online<double>(v.size()),
new moving_geometric_average_online<double>(v.size()),
new moving_harmonic_average_online<double>(v.size()),
new moving_generalized_average_online<double>(v.size()),
#else
new moving_arbitrary_generalized_average_online<double>(v.size(), identity_function<double>(), identity_function<double>()),
new moving_arbitrary_generalized_average_online<double>(v.size(), log, exp),
new moving_arbitrary_generalized_average_online<double>(v.size(), reciprocal_function<double>(), reciprocal_function<double>()),
new moving_arbitrary_generalized_average_online<double>(v.size(), square_function<double>(), sqrt),
#endif
};
double moving_results[] = {
#if !defined(USE_ARBITRARY_GENERALIZED_AVERAGE)
arithmetic_average(m.begin(), m.end()),
geometric_average(m.begin(), m.end()),
harmonic_average(m.begin(), m.end()),
generalized_average(m.begin(), m.end()),
#else
generalized_average(m.begin(), m.end(), identity_function<double>(), identity_function<double>()),
generalized_average(m.begin(), m.end(), log, exp),
generalized_average(m.begin(), m.end(), reciprocal_function<double>(), reciprocal_function<double>()),
generalized_average(m.begin(), m.end(), square_function<double>(), sqrt),
#endif
このように、関数ポインタや関数オブジェクトを指定することで、任意の一般化平均を求めることができるようになっています。但し、averages_online.hpp の方は C++11 または boost::function が必須となります。
さらに、statistics.hpp には、分散、標準偏差、絶対偏差、歪度、尖度などの各種統計量を、任意のコンテナ(配列、vector, list 等々)で以下のように求めるアルゴリズムが実装されています。
unbiased_variance(v.begin(), v.end()),
unbiased_standard_deviation(v.begin(), v.end()),
sample_variance(v.begin(), v.end()),
sample_standard_deviation(v.begin(), v.end()),
absolute_deviation(v.begin(), v.end()),
skewness(v.begin(), v.end()),
kurtosis(v.begin(), v.end()),
lag_one_autocorrelation(v.begin(), v.end()),
これらのうちいくつかの統計量についてのオンライン版のアルゴリズムはこの文献によって与えられています。そこで使われている式変形について多少補足しておきます。
平均:
\[\mu=E[X]\]
2次のモーメント:
\[\mu_2=E\left[\left(X-\mu\right)^2\right]=E\left[X^2-2X\mu+\mu^2\right]=E[X^2]-2E[X]\mu+\mu^2=E[X^2]-\mu^2=E[X^2]-E[X]^2=\sigma^2\]
3次のモーメント:
\[\mu_3=E\left[\left(X-\mu\right)^3\right]=E\left[X^3-3X^2\mu+3X\mu^2-\mu^3\right]=E[X^3]-3\mu(E[X^2]-E[X]^2)-\mu^3=E[X^3]-3\mu\mu_2-\mu^3\]
歪度:
\[\gamma=E\left[\left(\dfrac{X-\mu}{\sigma}\right)^3\right]=E\left[\dfrac{(X-\mu)^3}{\sigma^3}\right]=\dfrac{E[(X-\mu)^3]}{\sigma^3}=\dfrac{\mu_3}{\sigma^3}=\dfrac{E[X^3]-3\mu\mu_2-\mu^3}{\mu_2^{3/2}}\]
4次のモーメント:
\[\mu_4=E\left[\left(X-\mu\right)^4\right]=E\left[X^4-4X^3\mu+6X^2\mu^2-4X\mu^3+\mu^4\right]=E[X^4]-4E[X^3]\mu+6E[X^2]\mu^2-3\mu^4=\cdots=E[X^4]-4\mu\mu_3-6\mu^2\mu_2-\mu^4%\]
尖度:
\[\kappa=\dfrac{\mu_4}{\mu_2^2}-3=\dfrac{E\left[\left(X-\mu\right)^4\right]}{{\sigma^2}^2}-3=\dfrac{E[X^4]-4E[X^3]\mu+6E[X^2]\mu^2-3\mu^4}{\sigma^4}-3=\dfrac{E[X^4]-4\mu\mu_3-6\mu^2\mu_2-\mu^4}{\sigma^4}-3\]
さらにこれらの式により、移動平均と同様に移動分散、移動歪度、移動尖度などを求めることが出来ます。
それが実装された statistics_online.hpp の使用例として、その tarball にあるコードを以下に示しておきます。
/*
statistics.cc
Written by Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <vector>
#include <list>
#include <cstdlib> // random(), etc
#include <cassert>
#include "statistics.hpp"
#include "statistics_online.hpp"
int main()
{
const double maximum_permissible_error(1e-12);
const char *result_names[] = {
"mean",
"variance",
"standard deviation",
"biased variance",
"biased standard deviation",
"skewness",
"kurtosis",
"maximum value",
"minimum value",
"geometric mean",
"harmonic mean",
"root mean square",
};
std::vector<double> v(1000);
statistics_online<double> statistics_online;
srandomdev();
for (size_t i=0; i<v.size(); ++i) {
double x((double)random()/((1L<<31)-1));
v[i] = x;
statistics_online.add(x);
}
double results[] = {
arithmetic_average(v.begin(), v.end()),
unbiased_variance(v.begin(), v.end()),
unbiased_standard_deviation(v.begin(), v.end()),
sample_variance(v.begin(), v.end()),
sample_standard_deviation(v.begin(), v.end()),
skewness(v.begin(), v.end()),
kurtosis(v.begin(), v.end()),
*max_element(v.begin(), v.end()),
*min_element(v.begin(), v.end()),
geometric_average(v.begin(), v.end()),
harmonic_average(v.begin(), v.end()),
generalized_average(v.begin(), v.end(), 2),
};
double online_results[] = {
statistics_online.get_mean(),
statistics_online.get_unbiased_variance(),
statistics_online.get_unbiased_standard_deviation(),
statistics_online.get_sample_variance(),
statistics_online.get_sample_standard_deviation(),
statistics_online.get_skewness(),
statistics_online.get_kurtosis(),
statistics_online.get_maximum_value(),
statistics_online.get_minimum_value(),
statistics_online.get_geometric_mean(),
statistics_online.get_harmonic_mean(),
statistics_online.get_root_mean_square(),
};
for (size_t j=0; j<sizeof(online_results)/sizeof(online_results[0]); ++j) {
double error(sqrt(pow(online_results[j] - results[j], 2)));
std::cout << result_names[j] << ":\t" << online_results[j] << "\t" << results[j] << "\t" << error << std::endl;
assert(error < maximum_permissible_error);
}
const char *moving_result_names[] = {
"moving arithmetic average",
"moving geometric average",
"moving harmonic average",
"moving root mean square",
"moving sample variance",
"moving sample standard deviation",
"moving unbiased variance",
"moving unbiased standard deviation",
"moving skewness",
"moving kurtosis",
"moving minimum value",
"moving maximum value",
};
std::vector<double> m(v.size());
moving_statistics_online<double> moving_statistics_online(v.size());
for (size_t i=0; i<v.size()*10; ++i) {
double x((double)random()/((1L<<31)-1));
m[i % v.size()] = x;
moving_statistics_online.add(x);
if (i < v.size()) continue;
double moving_results[] = {
arithmetic_average(m.begin(), m.end()),
geometric_average(m.begin(), m.end()),
harmonic_average(m.begin(), m.end()),
generalized_average(m.begin(), m.end(), 2),
sample_variance(m.begin(), m.end()),
sample_standard_deviation(m.begin(), m.end()),
unbiased_variance(m.begin(), m.end()),
unbiased_standard_deviation(m.begin(), m.end()),
skewness(m.begin(), m.end()),
kurtosis(m.begin(), m.end()),
*min_element(m.begin(), m.end()),
*max_element(m.begin(), m.end()),
};
double online_moving_results[] = {
moving_statistics_online.get_mean(),
moving_statistics_online.get_geometric_mean(),
moving_statistics_online.get_harmonic_mean(),
moving_statistics_online.get_root_mean_square(),
moving_statistics_online.get_sample_variance(),
moving_statistics_online.get_sample_standard_deviation(),
moving_statistics_online.get_unbiased_variance(),
moving_statistics_online.get_unbiased_standard_deviation(),
moving_statistics_online.get_skewness(),
moving_statistics_online.get_kurtosis(),
moving_statistics_online.get_minimum_value(),
moving_statistics_online.get_maximum_value(),
};
for (size_t j=0; j<sizeof(online_moving_results)/sizeof(online_moving_results[0]); ++j) {
double error(sqrt(pow(online_moving_results[j] - moving_results[j], 2)));
if (i+1 == v.size()*10)
std::cout << moving_result_names[j] << ":\t" << online_moving_results[j] << "\t" << moving_results[j] << "\t" << error << std::endl;
assert(error < maximum_permissible_error);
}
}
return 0;
}
これら非オンラインとオンラインとの自乗誤差は 1,000 データにおいて 10E-12 程度に収まっていることは確認できています。
加えて、averages_online-20160119.tar.bz2 にて quantiles.hpp, quantiles_online.hpp も実装しました。これはメディアン等を求めるテンプレートライブラリです。しかし、オンライン版は移動平均しか現時点では実装しておりませんので、中途ゆえ解説は省略させて頂きます。使い方はソースコードを見れば一目瞭然だと思います。そして、averages_online-20160129.tar.bz2 にてヒューリスティック手法ですがオンライン版を実装しました。元になるアルゴリズムは、例えばこの文献などから辿れます。
それが実装された quantiles_online.hpp の使用例として、その tarball にあるコードを以下に示しておきます。
/*
quantiles.cc
Written by Taiji Yamada <taiji@aihara.co.jp>
*/
#include <iostream>
#include <vector>
#include <list>
#include <functional>
#if __cplusplus < 201103L
#include <boost/function.hpp>
#endif
#include <cstdlib> // random(), etc
#include <cassert>
#if __cplusplus < 201103L
#include <boost/utility.hpp>
#endif
#include "quantiles.hpp"
#include "quantiles_online.hpp"
int main()
{
const double maximum_permissible_error(1e-12);
const char *result_names[] = {
"p=0.25 quantile",
"p=0.5 quantile",
"p=0.75 quantile",
"0-quantile (least)",
"2-quantile (median)",
"1-quantile (greatest)",
};
std::vector<double> v(1000);
double probabilities[] = { .25, .5, .75, };
heuristic_quantiles_online<double> quantiles_online(probabilities);
srandomdev();
for (size_t i=0; i<v.size(); ++i) {
double x((double)random()/((1L<<31)-1));
v[i] = x;
quantiles_online.add(x);
}
sort(v.begin(), v.end());
double results[] = {
quantile_in_sorted_container(v.begin(), v.end(), 0.25),
quantile_in_sorted_container(v.begin(), v.end(), 0.5),
quantile_in_sorted_container(v.begin(), v.end(), 0.75),
least_in_sorted_container(v.begin(), v.end()),
median_in_sorted_container(v.begin(), v.end()),
greatest_in_sorted_container(v.begin(), v.end()),
};
double online_results[] = {
quantiles_online.get_quantile(0.25),
quantiles_online.get_quantile(0.5),
quantiles_online.get_quantile(0.75),
quantiles_online.get_minimum_value(),
quantiles_online.get_median(),
quantiles_online.get_maximum_value(),
};
for (size_t j=0; j<sizeof(online_results)/sizeof(online_results[0]); ++j) {
double error(sqrt(pow(online_results[j] - results[j], 2)));
std::cout << result_names[j] << ":\t" << online_results[j] << "\t" << results[j] << "\t" << error << std::endl;
if (!(j<=2 || j==4)) assert(error < maximum_permissible_error);
}
std::vector<double> m(v.size()), s(m.size());
moving_quantiles_online<double> moving_quantiles_online(v.size());
for (size_t i=0; i<v.size()*10; ++i) {
double x((double)random()/((1L<<31)-1));
m[i % v.size()] = x;
moving_quantiles_online.add(x);
if (i < v.size()) continue;
copy(m.begin(), m.end(), s.begin());
sort(s.begin(), s.end());
double moving_results[] = {
quantile_in_sorted_container(s.begin(), s.end(), 0.25),
quantile_in_sorted_container(s.begin(), s.end(), 0.5),
quantile_in_sorted_container(s.begin(), s.end(), 0.75),
least_in_sorted_container(s.begin(), s.end()),
median_in_sorted_container(s.begin(), s.end()),
greatest_in_sorted_container(s.begin(), s.end()),
};
double online_moving_results[] = {
moving_quantiles_online.get_quantile(0.25),
moving_quantiles_online.get_quantile(0.5),
moving_quantiles_online.get_quantile(0.75),
moving_quantiles_online.get_minimum_value(),
moving_quantiles_online.get_median(),
moving_quantiles_online.get_maximum_value(),
};
for (size_t j=0; j<sizeof(online_moving_results)/sizeof(online_moving_results[0]); ++j) {
double error(sqrt(pow(online_moving_results[j] - moving_results[j], 2)));
if (i+1 == v.size()*10)
std::cout << "moving " << result_names[j] << ":\t" << online_moving_results[j] << "\t" << moving_results[j] << "\t" << error << std::endl;
assert(error < maximum_permissible_error);
}
}
return 0;
}
さらに、averages_online-20160210.tar.bz2 にて、histogram.hpp, histogram_online.hpp を実装しました。
さらに、averages_online-20160302.tar.bz2 にて、multidimensional_histogram.hpp, multidimensional_histogram_online.hpp を実装しました。
さらに、averages_online-20160309.tar.bz2 の histogram.hpp, multidimensional_histogram.hpp にて、ヒストグラムからその確率密度に従う乱数を生成する histogram(_probability)_cumulative_distribution(_pseudo) を実装しました。
使い方は、ソースコードをご覧下さい。