C#からggplot2でグラフを描く
ggplot2は統計処理で用いられるR言語向けのグラフ作成ライブラリです。高機能かつ綺麗な出力が得られるので人気のようです。ただ、Rは使ったことがないので、今回はC#からこのライブラリを使ってみたいと思います。
R.NET
C#からですと当然Rをそのまま操作することはできませんが、R.Netを使うことで間接的にC#からRを使うことができます。
R.NET documentation -- user version | R.NET -- user version
R.NETの例
インストール方法はnugetを使うだけなので簡単です。事前にRを入れておく必要があります。
RからC#へ 整数の一様乱数を1000個生成
REngine.SetEnvironmentVariables(); using (REngine r = REngine.GetInstance()) { r.Initialize(); NumericVector x = r.Evaluate(@"runif(1000, min=1, max=10)").AsNumeric(); Console.WriteLine(x[0]); // 3.68815191206522 }
runifはRの乱数の生成関数で、Evaluateの引数がRのコマンドになります。変数xの配列には1000個の乱数が入っています。
C#からRへ 配列の合計を表示
using (REngine r = REngine.GetInstance()) { r.Initialize(); double[] value = { 1.1, 2.3, 3.5 }; NumericVector cv = r.CreateNumericVector(value); r.SetSymbol("rv", cv); // Rでの名前をつける double sum = r.Evaluate(@"sum(rv)").AsNumeric().First(); Console.WriteLine(sum); // 6.9 }
sumはRの合計関数です。
C#の変数をRの変数に変換することができればggplot2で表示できそうです。
ggplot2
このサイトを参考にして、散布図と近似曲線を表示してみます。
ggplot2 | R のグラフをより美しく
元のデータはC#で適当に作ります。横軸xが0から100の値で、縦軸yが横軸と同じ値に乱数を足したものです。
Rのxとyをそれぞれの軸に対応する変数名として定義しています。
var rd = new System.Random(); var x = Enumerable.Range(0, 200).Select(v => (double)v).ToArray(); var y = Enumerable.Range(0, 200).Select(v => (double)v + rd.Next(100)).ToArray(); // Rでの名前をつける r.SetSymbol("x", r.CreateNumericVector(x)); r.SetSymbol("y", r.CreateNumericVector(y));
ggplot2はデータをdataframe形式で受け付けるようなので変換が必要です。その後、データの指定、表示する形式を決めます。
ggplot2のオブジェクトgにどんどん表示したいものを足していくイメージのようです。
r.Evaluate("library(ggplot2)"); // dataframe形式 r.Evaluate("df <- data.frame(x=x, y=y)"); // データを指定 r.Evaluate(@"g <- ggplot(df, aes(x = x, y = y))"); // 表示形式を指定 // 散布図 r.Evaluate(@"g <- g + geom_point(size = 1.0)"); // 近似曲線の種類はmethodで選ぶ r.Evaluate(@"g <- g + geom_smooth(method = ""lm"")"); // ラベルの指定 r.Evaluate(@"g <- g + xlab(""x軸"")"); r.Evaluate(@"g <- g + ylab(""y軸"")"); r.Evaluate(@"g <- g + ggtitle(""タイトル"")"); // 表示 r.Evaluate(@"plot(g)");
劣モジュラ関数最小化のお勉強(Minimum Norm Point)
最近、劣モジュラ最適化と機械学習という本を買いました。その２章まで読んだので、劣モジュラ関数の最小化を試してみます。ちなみに２章は最小化、3章は最大化がテーマです。
劣モジュラ
集合を入力とする関数f()が劣モジュラ性(Sub modular)を持つのは以下の不等式を満たすときです。
$$
\begin{equation}
f(S)+f(T)\geq f(S \cup T)+f(S \cap T) (\forall S,T \subseteq V)
\end{equation}
$$
ここでVはとりうる最大の集合で、V={1}とか{1,2}とか{1,2,3}とかです。SとTはそれに含まれる集合なのでV={1,2,3}のときS={1,2}などがあり得ます。具体的な意味は[3]がわかりやすいです。イメージとしては大きな集合を入力したときほど、新たな要素を追加しても変化が小さい関数が劣モジュラ関数になります。逆に言うと小さい集合を入力しているときに、新たに要素を追加するとf()の値が大きく変わります。情報エントロピーもそのような関数で、あまり知らないときに受け取る情報の方が、よく知っているときより変化が大きいのです。
参考文献
- [1]劣モジュラ最適化と機械学習,2015.
- [2]SFO: A Toolbox for Submodular Function Optimization, 2010.
- [3]プリキュアで学ぶ劣モジュラ関数 - むしゃくしゃしてやった，今は反省している日記
- [4]The Minimum-Norm-Point Algorithm Applied to Submodular Function Minimization and Linear Programming, 2006.
- [5]Provable Submodular Minimization using Wolfe's Algorithm, 2014.
- [6]A note on the Minumum Norm Point algorithm, 2014.
最小化
劣モジュラ関数最小化は以下の式で表せます。
$$
\begin{equation}
\min f(S), S \subseteq V \tag{1}
\end{equation}
$$
関数fを最も小さくする集合Sを求めます。V={1,2,3}のときはS={},{1},{2},{3},{1,2},{1,3},{2,3},{1,2,3}の８通り(=$
\begin{equation}2^3=2^{|V|}\end{equation}$)あり、全部試せば良いわけです。ただしオーダーが指数関数なので、Vの集合か100種類とか1000種類とかになると相当時間がかかりそうです。
最小化にはタイトルにものせているMinimum-Norm-Point Algorithm(最小ノルム点アルゴリズム)を使うことができます。これを使うと全数探索より相当速く最小値の探索ができるようです。
L2ノルムの最小化
式(1)の最小化を行いたいのですが、最小ノルム点アルゴリズムではL2ノルムの2乗を最小化をすることでfの最小化を行います。
$$
\begin{equation}
\min \sum_{i=1}^n x_i^2,
x=(x_1,x_2,...,x_n) \in B(f) \tag{2}
\end{equation}
$$
ここでnはVのサイズで、xは後述するB(f)に含まれる点です。式(2)の最適解を$\begin{equation}x^*=(x_1^*,x_2^*,...,x_n^*)\end{equation}$とすると、fを最小化する集合は$\begin{equation}S^-= \big\{i \in V : x_i^*<0\big\} \end{equation}$となり[5]、式(2)を解くことは式(1)を解くことになります。
多面体
劣モジュラ関数fによって二つの多面体、劣モジュラ多面体$\begin{equation}P(f) \in \mathbb{R}^2\end{equation}$と基多面体$\begin{equation}B(f) \in \mathbb{R}^2\end{equation}$が定義できます。これは式(2)と関わりがあるので、最小化の計算に必要なものです。
劣モジュラ多面体
$\begin{equation}S \in V\end{equation}$について、$\begin{equation}x(S)\leq f(S)\end{equation}$をみたす$\begin{equation}x \in \mathbb{R}^2\end{equation}$の集合の全体が劣モジュラ多面体になります。ここで$\begin{equation}x(S)=\sum_{i \in S}x_i\end{equation}$とします。
$$
\begin{equation}
P(f)=\left \{ x \in \mathbb{R}^2:x(S)\leq f(S) (\forall S\subseteq V) \right \}
\end{equation}
$$
基多面体
劣モジュラ多面体に含まれるxのなかでも以下の式を満たす集合を基多面体と呼びます。
$$
\begin{equation}
B(f)=\left \{ x \in P(f):x(V)= f(V)\right \}
\end{equation}
$$
例
n=2、V={1,2}、S={1}or{2}or{1,2}のときは定義より以下の多面体になります。
$$
\begin{align}
P(f)=\left \{ x \in \mathbb{R}^2:x_1 \leq f(\{1\}),x_2 \leq f(\{2\}),x_1+x_2 \leq f(\{1,2\}) \right \} \\\
B(f)=\left \{ x \in \mathbb{R}^2:x_1 \leq f(\{1\}),x_2 \leq f(\{2\}),x_1+x_2 = f(\{1,2\}) \right \}
\end{align}
$$
f({1})=4,f({2})=3,f({1,2})=5を図示するとこんな感じで、B(f)は端っこの領域です。
端点集合
B(f)の端点（上図の黒丸）に注目してみます。n=2のときは2個、n=3のときは6個端点があります。これらはV={1,2,...,n}の成分番号集合を並び替えることによって、命名することが来ます。V={1,2}のときは2!=2なので、(1,2)と(2,1)があります。
この命名の仕方を線形順序と呼び、$\begin{equation}L=\left ( i_1,...,i_n \right )\end{equation}$とします。Lが決まると、B(f)の端点$\begin{equation}x^L=\left ( x_1^L,...,x_n^L \right )\in \mathbb{R}^n\end{equation}$が一意に決まります。L=(2,1)なら$\begin{equation}i_1=2,i_2=1\end{equation}$です。ここまで来ると$\begin{equation}x^L\end{equation}$を求めたくなってきます。
端点の導出
端点はGreedy法を使って求めることができます。この方法は軸ごとに順番に最大値を求めていくものです。$\begin{equation}x^{(2,1)}\end{equation}$ですとまずはじめに$\begin{equation}x_2\end{equation}$方向に最大値を探索します。上の図を見ると原点スタートで、上方向に探索を始めます。そうすると(0,3)で端にぶつかります。その後はそこから$\begin{equation}x_1\end{equation}$方向に進んでいき(2,3)で端点$\begin{equation}x^{(2,1)}\end{equation}$に到着します。$\begin{equation}x_1+x_2= 5\end{equation}$だったので、$\begin{equation}x_1+3=5\end{equation}$で2が導出されます。同様にはじめに$\begin{equation}x_1\end{equation}$方向に探索を進めれば、$\begin{equation}x^{(1,2)}\end{equation}$が得られます。
この処理を図形的解釈ではなく、定式化すると以下のようになります。
- $\begin{equation}L=(i_1,...,i_n),L(j)=\left\{i_1,...,i_j\right\}(j=1,...,n),L(0)=\{\}\end{equation}$
- $\begin{equation}x^L\end{equation}$の第$\begin{equation}i_j\end{equation}$成分を$\begin{equation}x_{i_j}^L=f(L(j))-f(L(j-1))\end{equation}$とする
ここで、L(j)はLの最初のj個の要素からなる部分集合です。
$\begin{equation}x^{\left (1,2,3 \right )}=\left (f(\{1\}),f(\{1,2\})-f(\{1\}),f(\{1,2,3\})-f(\{1,2\})\right )\end{equation}$です。
線形最適化
$\begin{equation}x\in B(f)\end{equation}$と$\begin{equation}c=(c_1,...,c_n)\in \mathbb{R}^n\end{equation}$の内積
$$
\begin{align}
\max \left \langle c,x \right \rangle, \\\
x=(x_1,x_2,...,x_n) \in B(f)
\end{align}
$$
これはV={1,...,n}の線形順序$\begin{equation}L=\left\{i_1,...,i_n\right\}\end{equation}$が$\begin{equation}c_1\geq c_1\geq ... \geq c_n\end{equation}$を満たすように取り、Greedy法を適用することでxを求めることができるそうです。$\begin{equation}c_1\end{equation}$から探索を始めるので、重要なものほど先に探索するイメージです。
VectorXd GetExtremePoint(const std::vector<int>& V, const VectorXd& c, std::function<double(const std::vector<int>&)> f) { // index付きのsortむずい std::vector<int> index(V.size()); std::iota(index.begin(), index.end(), 0); std::sort(index.begin(), index.end(), [&c](int a, int b){return c[a] > c[b];}); double pre_f = 0.0; VectorXd x(V.size()); std::vector<int> S; for (int j = 0; j < index.size(); j++) { int i_j = index[j]; S.push_back(V[i_j]); x[i_j] = f(S) - pre_f; pre_f += x[i_j]; } return x; }
Minimum Norm Point Algorithm
準備が長くなってしまいましたが、劣モジュラ関数の最小化法として速いと評判のMinimum Norm Point Algorithmを実装します。
Minimum norm pointをx*とするとf(S)が最小値となるのは、$\begin{equation}S=\{i|x^*_i<0\}\end{equation}$のときなので、x*が分かれば良いのです[5]。x*は$\begin{equation}\min_{x \in B(f)}||x||\end{equation}$の解です。
1. B(f)の端点を一つ求めて$\begin{equation}x^0\end{equation}$、$\begin{equation}\chi=\{x^0\}\end{equation}$,$\begin{equation}\hat{x}=x^0\end{equation}$とする
2. $\begin{equation}<\hat{x},x>\end{equation}$を最小化する$\begin{equation}x\end{equation}$を求める
$\begin{equation}<\hat{x},x>=\sum_{i}\hat{x}x_i\end{equation}$の最小化は$\begin{equation}c=-\hat{x}\end{equation}$として、Greedy法を使うと出てきます。その後$\begin{equation}\chi=\chi \cup \{x\}\end{equation}$をします。ここで最小値に近づいたら終了します。
3. $\begin{equation}\chi\end{equation}$のアフィン包$\begin{equation}aff(\chi)\end{equation}$において、$\begin{equation}\left \langle y,y \right \rangle\end{equation}$を最小化するyを求める
$$
\begin{align}
aff(S)=\left \{ y|y=\sum_{z\in S} \alpha_z z,\sum _{z \in S}\alpha_{z}=1;\alpha_z\geq 0\right \}\\\
conv(S)=\left \{ y|y=\sum_{z\in S} \alpha_z z,\sum _{z \in S}\alpha_{z}=1;\alpha_z\in \mathbb{R}\right \}
\end{align}
$$
aff(S)はアフィン包(affine hull)、conv(S)は凸包(convex hull)です。二点あったら、二点間の線分がconvで二点間の直線がaffのイメージです。
normの解は行列で解析的に解くことができて、Bを$\begin{equation}n\times |\chi|\end{equation}$行列とすると、$\begin{equation}\alpha=(B^\intercal B){\bf 1}/{\bf 1}^\intercal(B^\intercal B){\bf 1}\end{equation}$、$\begin{equation}y=B\alpha\end{equation}$です[5]。
yが$\begin{equation}conv(\chi)\end{equation}$の相対的内点なら$\begin{equation}\hat{x}=y\end{equation}$として1へ。
4. $\begin{equation}\hat{x}\end{equation}$とyの線分と$\begin{equation}conv(\chi)\end{equation}$の共通部分でyに近い点を探す
$\begin{equation}x'\in conv(\chi')\end{equation}$となるような極小な$\begin{equation}\chi'\end{equation}$を求めて2へ。
実行例
実装例
VectorXd AffineMinimizer(const std::vector<VectorXd>& Chi, MatrixXd &B) { int n = Chi[0].size(); B = MatrixXd::Zero(n, Chi.size()); for (int i = 0; i < Chi.size(); i++) { B.col(i) = Chi[i]; } VectorXd ones = VectorXd::Ones(Chi.size()); MatrixXd BtBt = (B.transpose() * B).inverse(); VectorXd alpha = (BtBt * ones) / (ones.transpose() * BtBt * ones); // 逆行列の計算をWoodburyの公式で減らすこともできるそうです return alpha; } VectorXd MinimumNormPoint(const std::vector<int>& V, std::function<double(const std::vector<int>&)> f) { // Step 1 // Chose any point x in B(f) and put X={x0}, x_hat=x0 VectorXd c(V.size()); VectorXd x0 = GetExtremePoint(V, c, f); std::vector<double> lambda; std::vector<VectorXd> Chi; Chi.push_back(x0); lambda.push_back(1.0); VectorXd x_hat = x0; while (1) { // Step 2 Find a point x that minimizes <x_hat, x> c = -x_hat; VectorXd x = GetExtremePoint(V, c, f); // if <x_hat, x> = <x_hat, x_hat> if (std::abs(x_hat.dot(x-x_hat)) < 0.0000001) { VectorXd x_star = x_hat; return x_star; } Chi.push_back(x); lambda.push_back(1.0); // Step 3 // Find the minumum norm point y in the affine hull of points in B(f) MatrixXd B; VectorXd alpha = AffineMinimizer(Chi, B); VectorXd y = B * alpha;std::cout << y << std::endl; if ((alpha.array() >= 0).all()) { x_hat = y; continue; } // Step 4 [5]を参考にしてます double theta = DBL_MAX; for (int i = 0; i < alpha.size(); i++) { if (alpha[i] < 0) { double temp_theta = lambda[i] / (lambda[i] - alpha[i]); if (theta > temp_theta) theta = temp_theta; } } x = theta * y + (1.0 - theta) * x; for (int i = 0; i < lambda.size(); i++) { lambda[i] = theta * alpha[i] + (1.0 - theta) * lambda[i]; } // list使って消すのが正解かも std::vector<double> new_lambda; std::vector<VectorXd> new_Chi; for (int i = 0; i < Chi.size(); i++) { if (lambda[i] > 0) { new_lambda.push_back(lambda[i]); new_Chi.push_back(Chi[i]); } } lambda = new_lambda; Chi = new_Chi; } } int main() { std::vector<int> V(10); std::iota(V.begin(), V.end(), 1); VectorXd x_star = MinimumNormPoint(V, EvaluateIwataFunction<10>); std::vector<int> S_star; for (int i = 0; i < x_star.size(); i++) { if (x_star[i] < 0) S_star.push_back(i); } }
劣モジュラ関数
Iwata's Test Function[4]
$\begin{equation}f(S)=|S||V \setminus S|-\sum_{j \in S}(5j-2n),(S\subseteq V)\end{equation}$
template <int n> double EvaluateIwataFunction(const std::vector<int>& X) { double sum = 0.0; for (int i = 0; i < X.size(); i++) sum += (5.0*(X[i])-2.0*n); return X.size() * (n - X.size()) - sum; }
本の例[1]
$\begin{equation}f(S)=2|S|(5 - |S|)- 8|S\cap \{1,3\}|+16\left |S\cap \{2\}\right |, (S\in \{1,2,3\})\end{equation}$
double EvaluateTestFunction2(const std::vector<int>& X) { double n = X.size(); double a = 2.0 * n * (5 - n); a += -8.0*((std::find(X.begin(), X.end(), 1) != X.end())?1:0) - 8.0*((std::find(X.begin(), X.end(), 3) != X.end())?1:0); a += 16.0*((std::find(X.begin(), X.end(), 2) != X.end())?1:0); return a; }
他にも情報エントロピーも劣モジュラ関数らしく、Classification問題で使えるようです。３章にのってるのかな。
Eigenの値をVisual Studioのデバッグモードで見る
Visual Studioのデバッグモードを使うと、実行中の変数の内容を見ることができます。ただし一般的な型以外は対応していないため見ることができません。私の場合ですと、Eigenという行列演算ライブラリーをよく使っているのですが、デバッグ中に値を見ることができないので不便に思っていました。
変数の内容表示
Visual Studio 2010以前では独自形式の変数表示を行う場合、autoexp.datを書き換えることで、独自形式の変数表示に対応します。
Developer's Corner - Eigen
一方、Visual Studio 2012以降では.XMLの.natvisファイルを使う方法に変わりました。知らずにautoexp.datをいじっても変化がなくてはまってしまいました。
やり方
eigen / eigen / source / debug / msvc / eigen.natvis — Bitbucket
このファイルをC:\Program Files (x86)\Microsoft Visual Studio 11.0\Common7\Packages\Debugger\Visualizersに格納します。
実行例
サーバーなしでD3.jsを使う
D3.jsはグラフの表示などができるJavaScriptのライブラリです。今回はこれを使って散布図を表示させてみます。
はじめに
データの可視化をする方法として一番有名なものはエクセルです。しかし、エクセルだと操作できるグラフを作ったり、かっこいいグラフを作るのは結構しんどいです。ほかに思いつく方法はPythonとかRとかMATLABですが、普通のパソコンには入っていないので、グラフ作成ツールを作っても共有できません。でもブラウザならどのパソコンにも入ってるし、便利かもと気づいたので、ブラウザで動くD3.jsを試してみました。
元ネタ
[1]» HTML5 の File API でドラッグ＆ドロップする TECHSCORE BLOG
[2]JavaScript - CSVをパースする - Qiita
[3]D3.js チュートリアル：基本5 – Axisで軸を表示して散布図を作る。 - Data is fun.
HTMLの準備
そもそもJavaScriptを使ったことがないので、勉強もかねています。（WEB系の技術も勉強しないと。。）
<html> <head> <title>ドラッグ＆ドロップで散布図</title> <script src="jquery-2.1.4.min.js"></script> <script src="d3/d3.min.js" charset="utf-8"></script> <script language="javascript">TODO</script> <body> <div id="droppable" style="border: gray solid 5px; padding: 2px; width:250px; height:50px;"> ここにファイルをドラッグ＆ドロップ </div> </body> </html>
このHTMLのdroppableなdivにファイルをドラッグ＆ドロップすると散布図を表示させるよう改造していきます。JavaScriptの制約でローカルにあるファイルにはアクセスできないですが、ドラッグ＆ドロップならオッケーみたいです。WEBサーバーを立てればローカルファイルでもいいんですけど、ソフトを配った人にサーバー入れて！なんて言えないですよね。
ドラッグ＆ドロップ
HTML5を使うとブラウザにファイルをドラッグ＆ドロップして読み込むことができます[1]。ファイルが読み終わると、drawScatter()という自作の関数を呼びます。この関数に散布図表示コードを書いていきます。
var fileReader = new FileReader(); fileReader.onload = function(event) { drawScatter(event.target.result); } fileReader.readAsText(file);
散布図の表示
データセット
表示させるデータはPRMLという本のclassificationを使います。前にもC#で試したことがありました。
wildpie.hatenablog.com
データ形式は1.208985 0.421448 0.000000みたいに、x軸 y軸 labelとなっています。labelは0と1の二種類です。これはcsvのパース関数を使うとそれぞれの数値を配列に入れることができます[2]。
function drawScatter(data) { var csvArray = parseCSV(data); console.log(csvArray); }
console.log()の結果を一部表示すると下記のようになります。これが200個ある感じです。
D3.jsで表示
D3.js チュートリアル[3]に散布図の作り方が詳しく載っていたので参考にしています。詳しくはこのサイトをご覧ください。
はじめに表示する領域を作ります。
var margin = {top: 20, right: 20, bottom: 30, left: 40}, width = 500 - margin.left - margin.right, height = 300 - margin.top - margin.bottom; var svg = d3.select("body").append("svg") .attr("width", width + margin.left + margin.right) .attr("height", height + margin.top + margin.bottom) .append("g") .attr("transform", "translate(" + margin.left + "," + margin.top + ")");
スケールの設定をします。データの範囲は[-2.5,3]でプロット画面は[0,width]の範囲とします。
var xScale = d3.scale.linear() .domain([-2.5, 3]) .range([0, width]); var yScale = d3.scale.linear() .domain([-3, 3]) .range([height, 0]);
軸の設定をします。
svg.append('g') .attr("class", "x axis") .attr("transform", "translate(0," + height + ")") .call(xAxis); svg.append('g') .attr("class", "y axis") .call(yAxis);
プロットするデータを定義します。csvArrayのデータを取り出して、'circle'を生成して、'cx'、'cy'で位置、'r'で大きさ、'fill'で色を指定します。
var colorCategoryScale = d3.scale.category10(); var points = svg.selectAll('circle') .data(csvArray) .enter() .append('circle') .attr('cx', function(d) { return xScale(d[0]) }) .attr('cy', function(d) { return yScale(d[1]) }) .attr('r', 7) .attr('fill', function(d) { return colorCategoryScale(d[2]); });
そんなこんなでこれを実行すると以下のようになります。
散布図を操作する
JavaScriptを使っているので、マウスで操作することができます。下の例ではマウスオーバー時は点を緑にして、大きさを２倍にしています。マウスオーバーをやめると元に戻ります。
points.on('mouseover', function() { d3.select(this) .attr('r', 14) .attr("fill", function(d) { return "green"; }); }) .on('mouseout', function() { d3.select(this) .attr('r', 7) .attr('fill', function(d) { return colorCategoryScale(d[2]); }) });
