C++で統計解析

C++ データマイニング 機械学習

先月開催されたJapan.R 2013の懇親会で，「C++で統計解析を行うための良いライブラリは？」という話がありました．

統計解析と一口に言っても結構広いので，ここでは以下の4つのカテゴリ

  • 記述統計量(最大値，最小値，平均値，分散等)
  • 統計的検定(t検定，χ2乗検定等)
  • 多変量解析(線形回帰，一般化線形モデル，判別分析，主成分分析，因子分析等)
  • 機械学習(サポートベクタマシン，ランダムフォレスト等)

に分けてライブラリがサポートする機能を整理してみると，下表のようになります*1．ここでは，Scytheなどの数値計算用のライブラリ，Shogun, Vowpal Wabbitなどの機械学習に特化したライブラリは対象外としています．他にもこんなライブラリがあるよ！という方は是非教えてください．

ライブラリ 記述統計量 統計的検定 多変量解析 機械学習
Boost.Accumulators
○
×
×
×
Apophenia
○
○
△
×
ROOT
△
△
△
△
GSL
△
×(?)
△
×(?)
ALGLIB
○
○
○
△

○：多くの手法を提供　△：少数の手法を提供　×：非提供

機能の網羅性という観点からはALGLIBが良さそうなので，本記事ではこのライブラリについて調べてみましょう．

ALGLIBとは

ALGLIBは，GPLライセンスのもとでOSSとして提供される数値計算およびデータ処理のためのソフトウェアです．C++, C#, Python, VBAなど多くのプログラミング言語から呼び出すことができます．商用版もあり，マルチスレッドIntel Math Kernel Libraryを用いた最適化・高速化などに対応しています．

以下のような機能を提供しています．

ダウンロードと解凍

以下のコマンドによりファイルをダウンロードし，解凍します．

$ wget http://www.alglib.net/translator/re/alglib-3.8.2.cpp.zip
$ unzip alglib-3.8.2.cpp.zip

cppというディレクトリが作成されます．

$ ls
cpp

ALGLIBが提供する機能は，cpp/srcディレクトリに格納されているヘッダファイルとcppファイルに記述されています．

$ ls cpp/src/
alglibinternal.cpp  dataanalysis.h      interpolation.cpp  solvers.h
alglibinternal.h    diffequations.cpp   interpolation.h    specialfunctions.cpp
alglibmisc.cpp      diffequations.h     linalg.cpp         specialfunctions.h
alglibmisc.h        fasttransforms.cpp  linalg.h           statistics.cpp
ap.cpp              fasttransforms.h    optimization.cpp   statistics.h
ap.h                integration.cpp     optimization.h     stdafx.h
dataanalysis.cpp    integration.h       solvers.cpp

データ構造

ALGLIBでは，以下のようなデータ構造が提供されています．

  • 扱える型
    論理値(boolean)，整数(integer)，倍精度浮動小数点(real)，複素数(complex)の4つ．複素数については，ALGLIB独自の実装も提供している．
  • 1次元の配列
    扱える型に応じて，boolean_1d_array, integer_1d_array, real_1d_array, complex_1d_arrayの4つ．
  • 2次元の配列
    扱える型に応じて，boolean_2d_array, integer_2d_array, real_2d_array, complex_2d_arrayの4つ．

使用例1) 相関係数の算出

まずはサンプルを参考に， x = (0, 1, 4, 9, 16, 25, 36, 49, 64, 81)^{\mathrm{T}}と y = (0, 1, 2, 3, 4, 5, 6, 7, 8, 9)^{\mathrm{T}}の間の相関係数を算出してみましょう．

#include <iostream>
#include "statistics.h"

using namespace std;
using namespace alglib;

int main()
{
  // ALGLIBを用いた相関係数の算出

  real_1d_array x = "[0,1,4,9,16,25,36,49,64,81]";
  real_1d_array y = "[0,1,2,3,4,5,6,7,8,9]";
  double v;

  cout.precision(8);

  // 共分散
  v = cov2(x, y);
  cout << "共分散: " << double(v) << endl;

  // Pearsonの相関係数
  v = pearsoncorr2(x, y);
  cout << "Pearsonの相関係数: " << double(v) << endl;

  // Spearmanの相関係数
  v = spearmancorr2(x, y);
  cout << "Spearmanの相関係数: " << double(v) << endl;

  return 0;
}

real_1d_arrayは1次元の倍精度浮動小数点の配列です．このデータ構造に，カンマ区切りの値を"[]"で囲んで文字列として代入しています．

srcディレクトリ内のすべてのヘッダファイルとcppファイル，および上記のソースコードのファイル(ここでは"test_cor.cpp"とします)を同じディレクトリに格納し，以下を実行します．

$ g++ -o test_cor *.cpp
$ ./test_cor 
共分散: 82.5
Pearsonの相関係数: 0.96269074
Spearmanの相関係数: 1

使用例2) 線形回帰

線形回帰は以下のように実行します(参照ページ)．

#include <iostream>
#include <vector>
#include "dataanalysis.h"

using namespace alglib;
using namespace std;

int main()
{
  // ALGLIBを用いた線形回帰

  // データセットの作成
  vector<double> xrange;
  vector<double> yrange;

  xrange.push_back(1.0); xrange.push_back(2.3); xrange.push_back(3.1);
  xrange.push_back(4.8); xrange.push_back(5.6); xrange.push_back(6.3);

  yrange.push_back(2.6); yrange.push_back(2.8); yrange.push_back(3.1);
  yrange.push_back(4.7); yrange.push_back(5.1); yrange.push_back(5.3);

  if (!xrange.empty() && !yrange.empty()) {
    // real_2d_arrayへのデータの格納
    int npoints = min(xrange.size(), yrange.size());//points
    real_2d_array xy;
    xy.setlength(npoints, npoints);
    for (int i = 0; i < npoints; i++) {
      xy(i, 0) = xrange[i];
      xy(i, 1) = yrange[i];
    }

    // 線形回帰の実行
    ae_int_t nvars = 1;
    ae_int_t info = 0;
    linearmodel lm;
    lrreport ar;
    lrbuild(xy,
            npoints,
            nvars,
            info,
            lm,
            ar);

    real_1d_array v;
    lrunpack(lm, v, nvars);

    // 回帰係数の出力
    // y = mx + b:
    // m (slope)  :  0.5842
    // b (y - int)  :  1.6842

    cout.precision(8);
    if (nvars > 0) {
      cout << "m: " << v(0) << " b: " << v(1) << endl;
    }
  }

  return 0;
}

real_2d_arrayにデータを格納して，lrbuild関数で線形回帰を実行し，lrunpack関数で回帰結果を取得しています．

コンパイルして実行します．

$ g++ -o test_lsfit *.cpp
$ ./test_lsfit 
m: 0.58418428 b: 1.6842239

使用例3) ランダムフォレスト

ALGLIBでは，Breiman(2001)のランダムフォレストのアルゴリズムを発展させたランダム決定フォレスト(Random Decision Forest)が実装されています．
インタフェースを調査中．分かり次第，書きます！

所感

統計解析は普段R(やPython)で行っているため，C++でこうした処理を行うこと自体が不思議な感じですが，やはり探索的なデータ解析には向いていないという印象です．
また，Rではデータフレームという型を中心として各種の統計解析の手法が実装されていますが，C++では特に複数のソフトウェアを用いる場合，共通した型がないためシームレスな連携が難しそうです．例えば，ALGLIBで主成分分析した結果をVowpal Wabbitで機械学習したいなどというニーズに応えるためには，型の変換や共通化が必要になるでしょう．

また機会があったら調べてみたいと思います．

*1:MCMCやベイズ的な手法など多くが漏れていますが，それらについてはまた今度別に整理したいところです