概要
ryouka0122さんの次の記事を見て閃きました。『これはSTLのいい練習問題になるのではないか』と……
数式とプログラミングの関係
※以下に登場する数式は上記記事より引用しています
※以下に使用する配列の変数名はdata、要素数はNとします
基本統計量
最大値・最小値 / max, min
\begin{align}
x_M&=\max_{i=1,...,N} {x_i}\\\
x_m&=\min_{i=1,...,N} {x_i}
\end{align}
#include <algorithm>
// max_elementは最大値のイテレータを取得する
// ゆえに*演算子で最大値の値に変換している
const auto max = *std::max_element(std::begin(data), std::end(data));
// 最小値を取得する
const auto min = *std::min_element(std::begin(data), std::end(data));
// 結果を表示する
std::cout << "最大値：" << max << std::endl;
std::cout << "最小値：" << min << std::endl;
総和 / Summation
S=\sum_{i=1}^{N}x_i
#include <numeric>
// accumulateは要素を加算した結果を返す
// (dataが浮動小数点数型なら第三引数を0.0にする必要がある)
const auto sum = std::accumulate(std::begin(data), std::end(data), 0);
std::cout << "総和：" << sum << std::endl;
平均 / Average, Mean
\bar{X}=\frac{1}{N}\sum_{i=1}^{N}x_i
#include <numeric>
// std::sizeはC++17以降で使用可能
// それ以前の場合、dataをstd::vectorやstd::arrayなどの
// コンテナで定義し、data.size()と書き換えること
const auto ave = std::accumulate(std::begin(data), std::end(data), 0.0) / std::size(data);
std::cout << "平均：" << ave << std::endl;
標本分散 / Variance
標本分散の元々の定義は1行目ですが、式変形により簡略化することができます。
しかし、一番下の式は「元のデータの二乗和」や「平均値の二乗」が大きな数値となった場合、計算誤差が元々の定義より大きくなる問題があります。Welfordの方法(漸化式を組む方法)によって解決できるそうですが、今回は素直に2通りのコードを書きます。
なお、「分散」には他に不偏分散もあり、「$\frac{1}{N}\sum_{i=1}^{N}(x_i-\bar{x})^2$」なところ、「$\frac{1}{N-1}\sum_{i=1}^{N}(x_i-\bar{x})^2$」にしたものとなります。なぜ2種類あるかの詳しい解説はこちら。
\begin{align}
Var&=\frac{1}{N}\sum_{i=1}^{N}(x_i-\bar{x})^2\\
&=\frac{1}{N}\sum_{i=1}^{N}(x_i^2-2x_i\bar{x}+\bar{x}^2)\\
&=\frac{1}{N}\sum_{i=1}^{N}x_i^2-\frac{2}{N}\sum_{i=1}^{N}x_i\bar{x}+\frac{1}{N}\sum_{i=1}^{N}\bar{x}^2\\
&=\frac{1}{N}\sum_{i=1}^{N}x_i^2-2\bar{x}\frac{1}{N}\sum_{i=1}^{N}x_i+\bar{x}^2\\
&=\frac{1}{N}\sum_{i=1}^{N}x_i^2-2\bar{x}^2+\bar{x}^2\\
&=\frac{1}{N}\sum_{i=1}^{N}x_i^2-\bar{x}^2\\
\end{align}
#include <numeric>
{
// 1つ目の手法(accumulate版) ※yumetodoさんのアイディア
// 第四引数に関数オブジェクトやラムダ式を渡せば加算以外の演算も可能
// (下コードでは引数のsumが直前の結果・eがその段階で読み込んだ配列の要素)
const auto ave = std::accumulate(std::begin(data), std::end(data), 0.0) / std::size(data);
const auto var = std::accumulate(std::begin(data), std::end(data), 0.0, [ave](double sum, const auto& e){
const auto temp = e - ave;
return sum + temp * temp;
}) / std::size(data);
std::cout << "分散：" << var << std::endl;
}
{
// 1つ目の手法(transform＆inner_product版)
// transformは各要素を処理・inner_productは内積計算
const auto ave = std::accumulate(std::begin(data), std::end(data), 0.0) / std::size(data);
std::transform(std::begin(data), std::end(data), std::begin(data), [ave](const auto &e){return e - ave;});
const auto var = std::inner_product(std::begin(data), std::end(data), std::begin(data), 0.0) / std::size(data);
}
{
// 2つ目の方法 ※akinomyogaさんのアイディア
const auto ave = std::accumulate(std::begin(data), std::end(data), 0.0) / std::size(data);
const auto var = std::inner_product(std::begin(data), std::end(data), std::begin(data), 0.0) / std::size(data) - ave * ave;
std::cout << "分散：" << var << std::endl;
}
{
// 2つ目の方法 ※_EnumHackさんのアイディア
// コンテナを1回走査するだけでOKにする
double ave = 0.0, var = 0.0;
for(const auto &x : data){
ave += x;
var += x * x;
}
ave /= std::size(data);
var = var / std::size(data) - ave * ave;
std::cout << "分散：" << var << std::endl;
}
標準偏差 / Standard Deviation
母標準偏差は標本分散の平方根、標本標準偏差は不偏分散の平方根です。式・コード略。
中央値 / Median
要するにデータ集団の真ん中を取ってくればいいだけですが、データ数が偶数なら隣接2項の平均値を取ります。
#include <algorithm>
std::sort(std::begin(data), std::end(data));
size_t median_index = std::size(data) / 2;
double median = (std::size(data) % 2 == 0
? static_cast<double>(data[median_index] + data[median_index - 1]) / 2
: data[median_index]);
std::cout << "中央値：" << median << std::endl;
最頻値 / Mode
値の範囲が決まっている場合は、別にプールを用意すれば$O(N)$の計算量で済みます。
※この項では便宜上、集計対象をunsigned intと仮定しています。
自動で判定させる場合、typename decltype(data)::value_typeとすれば集計対象の型が、
typename decltype(data)::iteratorとすれば集計対象のイテレーターが取得できます)
#include <algorithm>
// 集計する
std::vector<size_t> count(256, 0); //一例
for(const auto &x : data){
++count[x];
}
// 最大値の要素のインデックスを取り出す
auto max_iterator = std::max_element(count.begin(), count.end());
size_t mode = std::distance(count.begin(), max_iterator);
std::cout << "最頻値：" << mode << std::endl;
値が決まっていない場合は、std::unordered_map<値,カウント>としてからカウントの最大値を全走査するのが手っ取り早いでしょう。こちらも計算量は$O(N)$です(ハッシュが劣化してない場合)。
#include <algorithm>
#include <unordered_map>
// 集計する
std::unordered_map<unsigned int, size_t> hash;
for(const auto &x : data){
if(hash.find(x) != hash.end()){
++hash.at(x);
}else{
hash[x] = 1;
}
}
// 最大値の要素のインデックスを取り出す
// 別途比較関数を書きたくなかったのでラムダ式にした
// (ラムダ式の引数でautoが使えるのはC++14から)
auto max_iterator2 = std::max_element(hash.begin(), hash.end(),
[](const auto &a, const auto &b) -> bool {
return (a.second < b.second);
}
);
int mode = max_iterator2->first;
std::cout << "最頻値：" << mode << std::endl;
また、データを事前にソートして、先頭から見ていった際に「連続する数」をカウントして最頻値を探す作戦もあります(_EnumHackさんのアイディア)。
こちらはソートする時点で$O(N\log N)$になってしまいますが(カウント部分も結局は$O(N\log N)$になる)、「別途プールを用意する必要がない」といったメリットがあります。
※_EnumHackさんからのコメント：
はじめに、
std::adjacent_findで隣接する要素が等しい場所までイテレータを進めます。
次にその要素より値が大きい要素をstd::upper_boundで二分探索します。
最後にstd::distanceを使ってイテレータの距離（すなわち連続する数のカウント）を取得します。
これをコンテナの終端まで繰り返して最もカウントが大きかった要素をmodeとするという実装になっています。
別にデータ構造を用意（そしてそれにアクセス）するコストは馬鹿にならないと考えました。
そこで、あらかじめソートすることにより連続していない要素をstd::adjacent_findでつぎつぎに読み飛ばしながら、連続した要素を二分探索でカウントしつつイテレータを進めることで高速化を図る戦略を思いつきました。
#include <algorithm>
#include <vector>
#include <iostream>
std::sort(std::begin(data),std::end(data));
typename decltype(data)::value_type mode{};
size_t n{},count{1};
for(auto iter = std::adjacent_find(std::begin(data), std::end(data)),
last = std::end(data),
next = std::end(data);
iter != last;
){
next = std::upper_bound(iter, last, *iter);
count = std::distance(iter,next);
if( n < count ) n = count, mode = *iter;
iter = std::adjacent_find(next, last);
}
std::cout << "最頻値：" << mode << std::endl;
正規化 / Normalization
統計の分野における正規化(標準化)とは、ざっくり言えば次の式をデータ集団全体に施すことです(σは母標準偏差)。この処理によって、データ集団が「平均値0・標本分散1」となります。詳しい意味は次の記事などを参照ください。
統計学における標準化 - バイオインフォマティクス入門
\tilde{x}_i=\frac{x_i-\bar{x}}{\sigma}
#include <algorithm>
#include <vector>
// 既に平均値aveと標準偏差sdが計算されているとする
std::vector<double> norm_data(std::size(data));
std::transform(std::begin(data), std::end(data), norm_data.begin(), [&ave, &sd](const auto &e){
return (e - ave) / sd;
});
最後に
STLを使えばこんなにもコードがトリッｋ……ｹﾞﾌﾝｹﾞﾌﾝ、綺麗になることが分かりました。