C++で頻度分布を計算する

はじめに

C++でデータの頻度分布を計算するコードを紹介します。

ここではdoubleの変数をbinningして計算します。
C++特有のはまりどころや、log-binningするときのはまりどころもあるのでメモを残しておきます。

手順

方針は、doubleの値をbinningし、結果をmapに保持します。結果を出力するときは各binに入ったデータ点数をbinのサイズで割って出力します。
linear-scale, log-scaleの２種類を紹介しますが、基本的には同じ手順です。

normal-scaleのbinningの場合

コードはこちら

normal_bin.cpp

#include <iostream>
#include <cmath>
#include <map>
#include <random>

class HistoNormalBin {
public:
  HistoNormalBin(double bin_size) : bin(bin_size) {};
  void Add(double v) {
    int key = val_to_binidx(v);
    if( histo.find(key) == histo.end() ) { histo[key] = 0; }
    histo[key]++;
  }
  std::map<double,double> Frequency() {
    std::map<double,double> result;
    int key_min = histo.begin()->first;
    int key_max = histo.rbegin()->first;
    for( int i = key_min; i <= key_max; i++ ) {
      double val = binidx_to_val( i );
      size_t n = ( histo.find(i) == histo.end() ) ? 0 : histo[i];
      double freq = n / binidx_to_binsize( i );
      result[val] = freq;
    }
    return result;
  }
private:
  const double bin;
  std::function<int(double)> val_to_binidx = [=](double v)->int {
    return static_cast<int>( std::floor(v/bin) );
  };
  std::function<double(int)> binidx_to_val = [=](int i)->double {
    return i * bin;
  };
  std::function<double(int)> binidx_to_binsize = [=](int i)->double {
    return bin;
  };
  std::map<int, size_t> histo;
};

int main() {
  std::mt19937 mt(1234);
  std::uniform_real_distribution<double> dist(-50.0, 50.0);

  HistoNormalBin h(5.0);

  for( size_t i=0; i < 40; i++) {
    double r = dist(mt);
    h.Add(r);
  }

  auto histo = h.Frequency();
  for( const auto& keyval : histo ) {
    std::cout << keyval.first << ' ' << keyval.second << std::endl;
  }
}

HistoNormalBin が頻度を計算するためのクラスで、main関数の中で使用例を示しています。
(-50.0,50.0)の範囲で一様乱数を40回発生させて、その頻度分布を表示しています。

実装する上で注意すべき点は以下の通りです。

• データ点からbinのindex、binのindex->代表点、binのindex->binのサイズをそれぞれ計算するための関数を定義。
  • val_to_binidx, binidx_to_val, binidx_to_binsize という部分
• bin_idxを計算するときに、負の値が来てもいいように floor を使う。
  • 単純にintにcastするだけだと、例えば-0.5が0にキャストされてしまう。すると0に入るデータ点数が増えてしまう。
• C++のmapで、値があるかどうかを判定するときは histo.find(key) == histo.end() という式で判定する。Add()やFrequency()の中で、keyがない時にそのkeyの値を０にセットしている。
• Frequency()で結果を出力するときには、keyの最小値、最大値をとってその中でループを回すようにしている。
  • 頻度0のbinに対して0の値を入れるため。

log-scaleのbinningの場合

log-scaleのbinningの場合のコードはこちら

log_binning.cpp

#include <iostream>
#include <cmath>
#include <map>
#include <random>

class HistoLogBin {
public:
  HistoLogBin(double bin_base = 2.0) : base(bin_base) {};
  void Add(double v) {
    int key = val_to_binidx(v);
    if( histo.find(key) == histo.end() ) { histo[key] = 0; }
    histo[key]++;
  }
  std::map<double,double> Frequency() {
    std::map<double,double> result;
    int key_min = histo.begin()->first;
    int key_max = histo.rbegin()->first;
    for( int i = key_min; i <= key_max; i++ ) {
      double val = binidx_to_val( i );
      size_t n = ( histo.find(i) == histo.end() ) ? 0 : histo[i];
      double freq = n / binidx_to_binsize( i );
      result[val] = freq;
    }
    return result;
  }
private:
  const double base = 2.0;
  std::function<int(double)> val_to_binidx = [=](double v)->int {
    return static_cast<int>( std::floor( log(v)/log(base) ) );
  };
  std::function<double(int)> binidx_to_val = [=](int i)->double {
    return pow(base, i);
  };
  std::function<double(int)> binidx_to_binsize = [=](int i)->double {
    return pow(base, i);
  };
  std::map<int, size_t> histo;
};

int main() {
  std::mt19937 mt(12345);
  std::uniform_real_distribution<double> dist(0.0, 100);

  HistoLogBin h;

  for( size_t i=0; i < 10000; i++) {
    double r = dist(mt);
    h.Add(r);
  }

  auto histo = h.Frequency();
  for( const auto& keyval : histo ) {
    std::cout << keyval.first << ' ' << keyval.second << std::endl;
  }
}

normal binningの場合とほとんど同じです。
違うのはval_to_binidx, binidx_to_val, binidx_to_binsizeの３つのlambdaのみです。

おわりに

ここの手法を応用すれば、もっと複雑なbinningも容易に実装できます。
例えば、１以上のデータ点はnormal scale, １未満のデータ点はlog scaleでそれぞれbinningする、といったことも３つのlambdaを書き換えるだけで実装できます。