表計算で試すクラスタ分析
1. はじめに
今回の記事は、表計算ソフトで機械学習の一種であるクラスタ分析のk-means法を体験してみようというものです。
2. クラスタ分析手法：k-means法とは？
詳しく知りたい方は専門書やk-means法などで検索してヒットする専門的なサイトをご覧ください。
ここでは大まかにk-means法を解説します。
まず、何かの事象について2項目データが取得できたとします。それをX-Yの2次元散布図で表示してみると次の図の左上のようになっていたとします。直感的に「2つの群になっている！」と感じると思います。これがクラスタ分けをする動機です。そこで、クラスタ分析の1つ、k-means法を実施してみるとします。
この図の矢印に沿ってk-means法の処理順序を説明していきます。
- まずはこのデータ群1点1点にランダムに分類します。今回は2種類に分けようとしているので、○と×の2つを適当に割り当てます。
- ○と×それぞれの重心を求めます。重心の計算はｘ方向の平均、ｙ方向の平均です。それぞれの重心は◇と□とします。
- でたらめに決めた点の重心なので、それぞれの重心はほとんど同じ位置にありますが、拡大して見ると微妙にずれています。
- それぞれの重心に対して、各点の距離を算出し、近い重心に○と×を付け直します。つまり、◇に近ければ○にし、□に近ければ×に割当しなおします。
- 付け直した点で再び、○と×それぞれの重心を求めます。
- その新しい重心に対して、再び、各点の距離を算出し、近い方のしるしを付け直します。
- 区別がつくまで5に戻ります。区別がつくと思った時点で終わります。
以上のように、重心を求めて、近いものを割り当てる、を繰り返すのがk-means法と思って良いでしょう。
3. 表計算で試そう
利用するのは次です。
- 表計算ソフト（libreoffice-calc や Microsoft Excelなど）
- 次のタブ区切り文字列
今回はlibreoffice-calcの画面を使って説明します。
3.1. 表計算でのk-means法の内容
次の文字列をコピーします。
0 =MOD(A1,2) =(B1+1)*0.5+NORM.INV(RAND()*0.8+0.1,0,0.1) =(B1+1)*0.5+NORM.INV(RAND()*0.8+0.1,0,0.1) =(RAND()>0.5)*1 =1-E1 =E1*$C1 =E1*$D1 =F1*$C1 =F1*$D1 =SUM(G:G)/SUM(E:E) =SUM(H:H)/SUM(E:E) =SUM(I:I)/SUM(F:F) =SUM(J:J)/SUM(F:F) =($C1-K1)^2+($D1-L1)^2 =($C1-M1)^2+($D1-N1)^2 =IF(O1<P1,1,0) =1-Q1 =Q1*$C1 =Q1*$D1 =R1*$C1 =R1*$D1 =SUM(S:S)/SUM(Q:Q) =SUM(T:T)/SUM(Q:Q) =SUM(U:U)/SUM(R:R) =SUM(V:V)/SUM(R:R) =($C1-W1)^2+($D1-X1)^2 =($C1-Y1)^2+($D1-Z1)^2 =IF(AA1<AB1,1,0)
3.2. 表計算でk-means行う手順
表計算ソフトを起動し、A1にカーソルを合わせ、コピーした文字列をペーストします。
次のようになります。
セルA1の数を増やします。
増やしたセル分だけ、B1からAC1を下にコピーします
コピーした模様
足りないと思った場合は、再度A列を増やします。
増やした場合は、同じようにBからAC列も下に増やします。
ここから、クラスタ分析手法の1つ k-means法の繰り返し処理を行います。RからAC列をコピーします。
それをAD列にペーストします。
同じように、RからAC列を次々と右側にコピーしていきます。これが k-means法の繰り返し処理と同じことを行っています。
これでk-means法によるクラスタの推定が行われました。
3.3. データの図示
この分析がどんな形で行われたのかを見るために、グラフで図示しましょう。
まずは、CとD列を選択してグラフをプロットします。
グラフの種類は散布図を選びます。
プロットするとこのような点群が見えます。
人間がこの散布図を見たら、このデータは2つの群に分けられることが判ると思います。
しかし、計算機が分けるにはk-means法を使う訳ですが具体的に何をしたのかを見ましょう。
3.4. 算出結果の確認
算出結果の一番右から5列目のどれか1つを抜き出します。
次に、算出結果の一番右から7列目のどれか1つを抜き出し先ほどの下にその値が代入されるようにします。
一番右から6列目を先ほどの右上に
一番右から4列目を右下に
つまり、一番右下にからの列番号を示すと次のように作ります。
| X | Y | |
|---|---|---|
| a | 5つ右 | 7つ右 |
| b | 4つ右 | 6つ右 |
この4つのセルの値をグラフにプロットします。
グラフをダブルクリックし、さらにプロットエリアをダブルクリックし編集可能な状態にし右クリックし、データ範囲の編集を選びます。
追加のボタンをクリックします。
まずは X値をクリックします。
抽出した左側の値を選択します。
同じようにY値は右側の値を選択します。
範囲設定が終わったらOKを押します。
グラフに2つの点が現れているはずです。
この点をダブルクリックして、右クリックし線を追加します。
この現れた線が、k-means法で算出された二つのクラスタを分ける線です。
4. 解説
冒頭でk-means法というのは、重心を求めて、近いものを割当直す 手順であるという説明をしました。計算セル右側にコピーしていく作業のはこの繰り返しに相当します。計算セルの中身は次のようになっています。
全体的に次のような構成です。
- データ作成
- 割り当て初期化
- 重心計算
- 距離計算と再割り当て
4.1. データ作成
まずは、A,B,C,D列から。
| A | B | C | D |
|---|---|---|---|
| 0 | =MOD(A1,2) | =(B1+1)*0.5+NORM.INV(RAND()*0.8+0.1,0,0.1) |
=(B1+1)*0.5+NORM.INV(RAND()*0.8+0.1,0,0.1) |
A列は、ゼロから順に番号を割り当てています。
その隣のB列は、正解データを生成しています。
今回は、2種類の群を用意しますので、A列が偶数ならば0、奇数ならば1としています。
C,D列は、それぞれ2次元平面におけるX,Y値となるサンプルデータを生成しています。
サンプルデータは、群の基準位置にノイズを加えて生成します。
基準位置は、正解0の群では(x,y)=(0.5,0.5)にプロットし、正解1の群では(x,y)=(1,1)にプロットします。
ノイズは、一様分布の乱数RANDをNORM.INVを通し、ガウス分布の乱数にしたものを利用しました。
このようにして作られたテストデータはプロットすると、群の基準位置を中心にボンヤリとしたものになります。
4.2. 割り当て初期化
次に、E,F列について。
| E | F |
|---|---|
| =(RAND()>0.5)*1 | =1-E1 |
k-meansの初回は、群をでたらめに決めるところから始まります。
今回は、0.0から1.0まで一様分布の乱数を使い、0.5で2つに分けました。
F列は、E列の値を反転したものを格納しています。この値は後の重心計算に利用します。
4.3. 重心計算
その次、G,H,I,J と K,L,M,N列について。
| G | H | I | J | K | L | M | N |
|---|---|---|---|---|---|---|---|
=E1*$C1 |
=E1*$D1 |
=F1*$C1 |
=F1*$D1 |
||||
=SUM(G:G)/SUM(E:E) |
=SUM(H:H)/SUM(E:E) |
=SUM(I:I)/SUM(F:F) |
=SUM(J:J)/SUM(F:F) |
G,H,I,J列は、X値の0群と1群、Y値の0群と1群の値を格納しています。それぞれ属してない群については0にしています。こうして0にすることで、例えば、G列全部の合計は、0群のX値の合計に等しくなります。
これを0群の個数で割れば、平均値を求められます。個数は、先ほどのE列の合計と等しいです。反転したものを算出しておいたのは、このためです。
この計算をK,L,M,N列で行っています。
なお、これらの値は列全部同じ値になるので、1つで十分なのですが、下方向のコピーだけで次の距離計算を行いたかったため、このように下にコピーしました｡
4.4. 距離計算と再割り当て
最後に、O,P,Q,R列について。
| O | P | Q | R |
|---|---|---|---|
=($C1-K1)^2+($D1-L1)^2 |
=($C1-M1)^2+($D1-N1)^2 |
=IF(O1<P1,1,0) |
=1-Q1 |
O,P列は、各点の重心からの距離を算出しています。中学校で習ういわゆるピタゴラスの定理、三平方の定理とほぼ同じなのですがちょっと違うのは最後に平方根をとってない点です。ここで距離を算出する目的は大小を比較するためだけですので、平方根をとらなくても結果は同じになります。
Q列では、その距離の低い方の重心に属する方に再割り当てしています。R列は、Q列の反転です。次の計算で、属する個数を計算しやすくするためです。
そしてこれ以降は、Q,R列を使って、再び重心計算し、再び距離計算し、再割り当てし、そして、再び重心計算し。。。と繰り返します。
今回の例では、2回繰り返すだけで、正しい群に分けられています。
5. (おまけ)ヒストグラムで見るクラスタ
2次元の散布図で見ると分かるクラスタですが、1次元のヒストグラムではどんな形になっているかを見てみましょう。
まずは空いている所に、0, 0.1と縦に作りその2つを選択し右下の黒いぽっちをドラッグして1.2まで作りましょう。
その隣りに同じ行数選択して数式を入力しましょう。
先ほどの数字をセルAE1からAE13に作った場合のものですので、違う場合はその部分を読み替えて下さい。始めと終わりの波かっこが重要です。
{=frecuency( c:c, ae1:ae13 )}
これを入力し終えたら、enterを押さずに、Ctrl+Shift+Enterで入力完了します。
すると、ヒストグラムが計算されて表示されます。
あとはこれらを棒グラフで表示するだけです。
棒グラフを選びましょう。
最初の列を項目名にすると意図したグラフが作られます。
ヒストグラムが表示されました。今回のは1次元でヒストグラムをとっても2つのクラスタが現れているのがよく分かります。
今回のデータは、ガウス分布のノイズを加えており、分布としては2つの山が見えるものになっています。
実際に何かのデータ取得して、ヒストグラムに2つの山が見えたら、何か別々の事象が起きていることを疑うのが良いと思います。目的のデータが正しく取得できているのか？測定対象に何か問題が起きてないか？ということを考えるきっかけになります。
なお、ノイズをガウス分布ではなく、一様分布にするとヒストグラムでは群に分けることが難しくなることがあります。この場合、多次元を一度に演算できるk-means法がとても有効になってきます。
6. おわりに
今回は、表計算ソフトを用いて機械学習の流れをみてみました。これを一通り試すと普通のプログラムで作った方が楽だと感じるのではないでしょうか？
しかし、それこそが理解に近づいたという証拠だろうと思います。
現代の人工知能は機械学習がベースになっています。理解するためにこのようなアプローチも一つだろうと思います。
Discussion