思考ノイズ

無い知恵を絞りだす。無理はしない。

<はじパタ> 10.2 非階層型クラスタリング (K-平均法)

いろいろごたごたしている中でも少しづつはじパタすすめておりまして、とりあえず10章のストックができたのでこちらの備忘録を掲載します。

10章はクラスタリングで教師データなしの状態でパターンニングをおこなう手法になります。定番のアヤメデータをこのK-平均法でクラスタリングをおこなうとどうなるかをスクリプトを書いてみてみました。

アヤメデータ

このアヤメデータは3種50セット、合計150個のアヤメそれぞれの特性(がく片の幅、がく片の長さ、弁の幅および 弁の長さ)を記載しています。

例えば、がく片の幅と弁の長さでデータをプロットした場合の2次元データは以下のプロットで表現できます。

f:id:bython-chogo:20180403003147p:plain

つまり、今回の問いはこの150個のデータを教師データなしでとりあえず3種区別してみる、というものになります。 アヤメデータのプロットのスクリプトは以下のサイトを参考にさせていただいています。

Python: Iris データセットを Matplotlib で可視化してみる | CUBE SUGAR STORAGE

今回のスクリプトについて

今回は視覚化しやすいように4つの特性から2つをとったデータについてK平均法を使って3つのクラスタに分けてみます。4つのデータから2つとるので、6種のデータとなります。さらに4つのデータ全部を使ってどう変化するかもみてみたいと思います。

評価方法

評価には2つの指標を使います。ひとつは正解率、もうひとつは精度と再現率を評価できるF値を求めます。以下のリンクが詳しいかと思います。

抑えておきたい評価指標「正解率」「精度」「再現率」 - HELLO CYBERNETICS

K平均法

K-平均法は最初にそれぞれの値のクラスタをランダムに決定します。今回は3個のアヤメデータなので3つのクラスタにわけます。そしてその3つのクラスタの平均値のポイントを求め、またそれぞれの値で一番距離が近い平均値のクラスタに変更し直します。これを繰り返してデータが更新されなくなったら終わりとするものです。

実行結果

がく片の幅と弁の長さの結果

まずは先ほども例として挙げた「がく片の幅」と「弁の長さ」についてK-平均法の結果をプロットしてしてみます。

  • K-平均法の結果

f:id:bython-chogo:20180403005841p:plain

  • 実際のデータ

f:id:bython-chogo:20180403131050p:plain

  • 評価値
 ===  sepal width (cm) ,  petal length (cm)  === 
CORRECTNESS:  0.3867 ,   F-Value    :  0.4761
CORRECTNESS:  0.94 ,     F-Value    :  0.7417
CORRECTNESS:  0.8467 ,   F-Value    :  0.8337
CORRECTNESS:  0.9 ,  F-Value    :  0.897
CORRECTNESS:  0.9333 ,   F-Value    :  0.9328
CORRECTNESS:  0.9533 ,   F-Value    :  0.954
CORRECTNESS:  0.96 ,     F-Value    :  0.9607
CORRECTNESS:  0.94 ,     F-Value    :  0.9413
CORRECTNESS:  0.9267 ,   F-Value    :  0.9271

もとのデータが比較的わかりやすく分かれていることもあるので、正答率もわるくないと思われます。ランダムにはじめて平均値とその距離でこうもきれいにわかれるものなのかと少し驚きました。3つの平均値をみてみると最初は3つとも近い場所にいるのですが、繰り返していくうちに平均的なポジションに落ち着く様子がわかりました。これも今度アニメーションでプロットしてみたいと思います。

また正解率やF値はピークを越えて少しだけ下がるところで収束してしまうようです。実際は回答とすべき教師データはないので、収束まつしか判断のしようがなさそうです。

4つの特性から2つピックアップしてみる

では、4つの特性から2つピックアップした時のそれぞれのK-平均法の評価値はどのように変わるか見てみましょう。初期値がランダムなので評価ごとに値は変更しますが、とりあえず評価値の大小関係は変わらないと考えてよいです。

 ===  sepal length (cm) ,  sepal width (cm)  === 
CORRECTNESS:  0.4 ,  F-Value    :  0.4831
CORRECTNESS:  0.7867 ,   F-Value    :  0.689
CORRECTNESS:  0.8067 ,   F-Value    :  0.8134
CORRECTNESS:  0.8267 ,   F-Value    :  0.8363
CORRECTNESS:  0.82 ,     F-Value    :  0.8306

 ===  sepal length (cm) ,  petal length (cm)  === 
CORRECTNESS:  0.3867 ,   F-Value    :  0.4652
CORRECTNESS:  0.8267 ,   F-Value    :  0.7216
CORRECTNESS:  0.88 ,     F-Value    :  0.8794
CORRECTNESS:  0.9 ,  F-Value    :  0.9021
CORRECTNESS:  0.8867 ,   F-Value    :  0.8918
CORRECTNESS:  0.88 ,     F-Value    :  0.885
CORRECTNESS:  0.8867 ,   F-Value    :  0.8897
CORRECTNESS:  0.88 ,     F-Value    :  0.8824

 ===  sepal length (cm) ,  petal width (cm)  === 
CORRECTNESS:  0.3667 ,   F-Value    :  0.4585
CORRECTNESS:  0.8333 ,   F-Value    :  0.6973
CORRECTNESS:  0.8267 ,   F-Value    :  0.83
CORRECTNESS:  0.8267 ,   F-Value    :  0.8352
CORRECTNESS:  0.8333 ,   F-Value    :  0.8414

 ===  sepal width (cm) ,  petal length (cm)  === 
CORRECTNESS:  0.3867 ,   F-Value    :  0.4761
CORRECTNESS:  0.94 ,     F-Value    :  0.7417
CORRECTNESS:  0.8467 ,   F-Value    :  0.8337
CORRECTNESS:  0.9 ,  F-Value    :  0.897
CORRECTNESS:  0.9333 ,   F-Value    :  0.9328
CORRECTNESS:  0.9533 ,   F-Value    :  0.954
CORRECTNESS:  0.96 ,     F-Value    :  0.9607
CORRECTNESS:  0.94 ,     F-Value    :  0.9413
CORRECTNESS:  0.9267 ,   F-Value    :  0.9271

 ===  sepal width (cm) ,  petal width (cm)  === 
CORRECTNESS:  0.3933 ,   F-Value    :  0.4796
CORRECTNESS:  0.9533 ,   F-Value    :  0.7596
CORRECTNESS:  0.86 ,     F-Value    :  0.7473
CORRECTNESS:  0.8467 ,   F-Value    :  0.841
CORRECTNESS:  0.8733 ,   F-Value    :  0.876
CORRECTNESS:  0.8933 ,   F-Value    :  0.8965
CORRECTNESS:  0.9133 ,   F-Value    :  0.9168
CORRECTNESS:  0.9267 ,   F-Value    :  0.9293
CORRECTNESS:  0.9267 ,   F-Value    :  0.9292

 ===  petal length (cm) ,  petal width (cm)  === 
CORRECTNESS:  0.36 ,     F-Value    :  0.4523
CORRECTNESS:  0.9467 ,   F-Value    :  0.7468
CORRECTNESS:  0.8467 ,   F-Value    :  0.8292
CORRECTNESS:  0.9267 ,   F-Value    :  0.9258
CORRECTNESS:  0.96 ,     F-Value    :  0.9602
CORRECTNESS:  0.9733 ,   F-Value    :  0.9735
CORRECTNESS:  0.9667 ,   F-Value    :  0.9672
CORRECTNESS:  0.96 ,     F-Value    :  0.9606

概ね、0.85~0.95といったあたりに収まるでしょうか。

4つの特性を使った評価

では、2つといわず一気に4つの特性を使った4次元での評価はどうなるか確認をしてみたいと思います。

CORRECTNESS:  0.4267 ,    F-Value    :  0.5037
CORRECTNESS:  0.8467 ,   F-Value    :  0.7439
CORRECTNESS:  0.8667 ,   F-Value    :  0.8586
CORRECTNESS:  0.9 ,  F-Value    :  0.897
CORRECTNESS:  0.9 ,  F-Value    :  0.9007
CORRECTNESS:  0.92 ,     F-Value    :  0.9221
CORRECTNESS:  0.92 ,     F-Value    :  0.9227
CORRECTNESS:  0.9067 ,   F-Value    :  0.9096
CORRECTNESS:  0.9 ,  F-Value    :  0.9021
CORRECTNESS:  0.8933 ,   F-Value    :  0.8937
CORRECTNESS:  0.8867 ,   F-Value    :  0.8866

あれれ、精度としては特徴がよく表れている2つをピックアップしたほうがよい数値がでているようにもみえます。でもおそらく、その結果はその特徴量に強く適合がしすぎてしまっているかと思われるので、こちらの4次元での方法はより「一般的」「総合的」な評価がしやすくなるのかとは思います。

ソースコード

今回利用したソースコードです。Jupyter notebook上で動作を確認しており、 processで実行、2つの引数はアヤメの4つ特性 0~3のうち2つピックするのと、3つ目の引数はプロットをするかTrue/Falseで表しています。

import numpy as np
from sklearn import datasets
import matplotlib.pyplot as plt
from collections import Counter
#inline pyplot

class Iris_PM:
    def __init__(self):
        self.p = [2, 3]
        self.x = []
        self.y = []
        self.t = []
        self.pt = []
        self.num = [50, 50, 50]
        iris = datasets.load_iris()
        features = iris.data
        targets = iris.target
    def install(self, a=2, b=3):
        self.x = features.T[a]
        self.y = features.T[b]
        self.t = targets

    def plot_map(self, df=False):
        p = self.pt
        if df:
            p = self.t

        for t, marker, c in zip(range(3), '>ox', 'rgb'):
             
            plt.scatter(
                #self.x[50*t:50*t+49],
                #self.y[50*t:50*t+49],
                features[p == t, self.p[0]],
                features[p == t, self.p[1]],
                marker=marker,
                c=c,
            )
            plt.xlim(1.9, 4.5)
            plt.ylim(0.9, 7.1)
            plt.xlabel(feature_names[self.p[0]])
            plt.ylabel(feature_names[self.p[1]])
            plt.autoscale()
            plt.grid()
        #print(linear_pro(w, np.array([1.,1.,1.])))
        plt.show()
    def randomize(self):
        self.pt = np.array([np.random.randint(3) for i in range(sum(self.num))])
    
    def means(self):
        nu = [[0.0, 0.0], [0.0, 0.0], [0.0, 0.0]]
        cnt = [0, 0, 0]
        for i in range(sum(self.num)):
            nu[krs.pt[i]][0] += self.x[i]
            nu[krs.pt[i]][1] += self.y[i]
            cnt[krs.pt[i]] += 1
        for i in range(3):
            nu[i][0] /= cnt[i]
            nu[i][1] /= cnt[i]

        return nu
    def correct_rate(self):
        #print(self.pt)
        cnt = np.array([[Counter(self.pt[i*50:i*50+50])[j] for i in range(3)] for j in range(3)])
        #print(cnt)
        #print([np.argmax([Counter(self.pt[j*50:j*50+50])[i] for i in range(3)]) for j in range(3)])
        #print([np.argmax([Counter(self.pt[i*50:i*50+50])[j] for i in range(3)]) for j in range(3)])
        odr = [np.argmax([Counter(self.pt[j*50:j*50+50])[i] for i in range(3)]) for j in range(3)]
        lodr = []
        for i in odr:
            lodr += ([i]*50)

        fv = 1.0
        
        for i in range(3):
      
            x = np.count_nonzero(self.pt[i*50:i*50+50] == odr[i])
            z = np.count_nonzero(self.pt[(i+1)%3*50:(i+1)%3*50+50] == odr[i]) + \
                np.count_nonzero(self.pt[(i+2)%3*50:(i+2)%3*50+50] == odr[i])
            rec = x/50
            pre = 50/(50+z)
            fv *= 2*rec*pre/(rec+pre)

        print("CORRECTNESS: ", "{:.4}".format(np.sum(self.pt == lodr)/150), ",\t F-Value    : ", "{:.4}".format(np.cbrt(fv)))
        
    def process(self, a=1, b=2, plot_map=False):
        self.install(a,b)
        self.randomize()

        bak_list = []
        cnt = 0
        print(" === ", feature_names[a], ", ", feature_names[b]," === ")
        while self.pt is not bak_list and cnt < 30:

            self.correct_rate()

            cnt += 1
            bak_list = self.pt
            nu = self.means()

            up = []
            #krs.plot_map()
            for i in range(sum(self.num)):
                #print(n,":", krs.x[i], krs.y[i], nu)
                up.append(np.argmin([[(self.x[i]-n[0])**2 + (self.y[i]-n[1])**2 for n in nu]]))

            self.pt = np.array(up)
            nu = self.means()


            if np.allclose(bak_list, self.pt):
                if plot_map:
                    self.plot_map()
                    self.plot_map(True)
                break        
        print("")
krs = Iris_PM()

krs.process(1,2,False)