思考ノイズ

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

反復性肩脱臼の手術から1年経過

早いもので反復性肩脱臼からおよそ一年たちました。

bython-chogo.hatenablog.com

このブログでは機械学習系の記事をメインで書いていますが、いまだに半分以上のアクセスがこちらの記事のようです。同じ悩みを抱える人は意外と多く、どのような手術なのか気になって人も多いみたいですね。こういうレポートって大事なんだなとアクセスログを見るたびに思います。そこで一年後どうなっているかの記録と、思い出せる限りの1年の経過も残しておきたいと思います。

現在、1年後の状況

当然ですが、不自由なく普通に生活ができています。手術した左肩はきちんと後ろまでまわせますし、脱臼することはないです。可動域は右と比べて少ないのは否めないと思います。後ろに回すと引っ張られる感じはありますが、人口のひもで固定しているので仕方ないのかと思います。あとリハビリを頑張ればもっと行くのかもしれませんが、私生活でその必要は感じられないので私個人としては十分といえます。

一年の経過

  • 術後すぐ:術部の回復をおこなうため、サポータによる固定で過ごします。利き手ではなかったので、まだよかったですが、それでも不自由でした。寝るときもつけたままでねるので、寝返りもできないし、ずっと吊っている左肩の肩こりはとてもひどかったのを覚えています。リハビリも週に2回です。
  • 術後一か月経過:サポーターは外せますが、可動域はとても狭いです。前から上に手を上げようとしても最初は90度まであがりませんでした。生活にも影響がでます。左手を使って胴体から離して行わなくてはいけない作業は大変です。高いものを取るとかは右手でもできますが、何気に靴ひもを結ぶことは両手を使わなくてはならず、さらに胴体から離れてやる作業なので難しいです。リハビリは引き続き週二回です。
  • 術後三か月経過:だんだんと可動域が広がっていきますがまだ制限がでます。前から上に腕を上げる運動は180度までもう少しのところで止まってしまい、なかなかそこから改善せずもどかしい思いをしました。診断の時もリハビリが芳しくないということでもし次回の診断までによくなければ筋肉を柔らかくする注射をうつといわれました。結果として一か月後の診断では合格点がでてその注射はうたなくてもよくなりましたが。それを受けてリハビリの回数も週1回になりました。
  • 術後六か月経過:運動が解禁されて、ジム通いもスタートしました。可動域も実生活に問題がない程度まで動くようになりました。そこから数か月してリハビリも終わりとなりました。

というわけでリハビリが不安になった時もありましたが、結果として日常生活にはなにも問題ないレベルまで回復しました。脱臼するかもしれないと可動域をセーブしていた術前が嘘のように特に何も気にせずに生活をしています。半年以上のリハビリはありましたが、今後の人生や老後のことも考えるとやってよかったとはっきりと言えます。お悩みの方は半年以上はすこし不自由ですが、なるべく早いうちに決断することをおすすめします。

f:id:bython-chogo:20180501002752p:plain:w300

<SVM> サポートベクタマシン自主勉強(1):方針

ちょっと間が空いてしまいましたが機械学習のエントリとなります。前回最後にちらっとお話ししたのですが、サポートベクタマシン(SVM)に足を突っ込んでみたいと思っています。

SVMのコーディングが自力でできない問題

ただ、このSVMは何回か挑戦して失敗している経緯があります。前回のようにはじパタ本だけでは私のレベルでクリアができない代物となります。そこではじパタと合わせて、「異常検知と変化検知」の内容と相互補完してみています。

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)

異常検知と変化検知 (機械学習プロフェッショナルシリーズ)

そこでやっと、SVMの概要について理解した。。。。つもりでした。しかしながらいざこれをスクリプトに落とそうとするとどう書いていいのかわからない。つまりロジックがわかったつもりでも、実行に移すとなるとどうしていいのかわからなくなってしまいました。

頼るべきはGoogle先生集合知

しかしながらもうすでに歴史のあるこの分野では検索してしまえばいくらでもコードが出てきます。今回は以下のサイトさまのコードを拝借してやはりアヤメデータを使った分類を行ってみました。

qiita.com

こちらのコードは y = x の線形で分離をして  t を1, -1とわけています。アヤメデータもとりあえず、3種のうちの最初の1種を1, その他を-1としてインプットするコードに改造をしてみました。 その結果が以下となります。

f:id:bython-chogo:20180430230944g:plain

おぉ、なんとなく分類する線が学習を経るごとに赤と緑の間にはいっていきますね。赤が1としたアヤメデータなのでうまくいっているようです。

ここからスタート

SVMについてまとめることはたくさんありそうです。ラグランジュ乗数法、正則化、双対問題、KKT条件、カーネルトリック。数学的要素が多くどこまでカバーできるかわかりませんが、まず次回はこのスクリプトについて解いていくことから初めてみようかと思います。正直きちんと着地できるかまだ不明瞭ではあります。

コード

Jupyter Notebookで実行しています。

import sys
import numpy as np
from sklearn import datasets
import itertools
import matplotlib.pyplot as plt
#inline pyplot

from matplotlib import animation
from IPython.display import HTML

iris = datasets.load_iris()
features = iris.data
feature_names = iris.feature_names
targets = iris.target

class SVM():
    def __init__(self, v1=0, v2=2, t=0):
        self.X = np.array([features.T[v1], features.T[v2]]).T
        self.T = np.array([1 if i==t else -1 for i in targets])
        self.N = self.X.shape[0]
        self.d = self.X.shape[1]
        self.alpha = np.zeros(self.N)
        self.beta = 1.0
        self.eta_al = 0.0001
        self.eta_be = 0.1
        self.v = [v1, v2]

    def test(self):
        print(self.X.shape)
        print(self.T)
        print(self.N)
        print(self.alpha)
        
    def plot_map(self, w, b):
        plt.subplot(1, 1, 1)
        seq = np.arange(4.0, 8.0, 0.02)

        plt.xlabel(feature_names[self.v[0]])
        plt.ylabel(feature_names[self.v[1]])
        plt.xlim(4, 8)
        plt.ylim(0, 8)

        plt.plot(seq, -(w[0] * seq + b) / w[1], 'k-')
        #plt.autoscale()
        plt.grid()
        for t, marker, c in zip(range(3), '>ox', 'rgb'):
            #print(self.X[self.T==1], 0)
            #print(self.T)
            plt.scatter(
                features[targets == t, self.v[0]],
                features[targets == t, self.v[1]],
                marker=marker,
                c=c,
            )
            #print(linear_pro(w, np.array([1.,1.,1.])))
        plt.show()
        
        
    def rtin(self, show=False):
        
        k_w = np.array([-100.0, -100.0])
        k_b = -100.0
        value = []
        for _itr in range(1000):
            for i in range(self.N):

                self.delta = 1 - (self.T[i] * self.X[i]).dot(self.alpha * self.T * self.X.T).sum() - self.beta * self.T[i] * self.alpha.dot(self.T)
                self.alpha[i] += self.eta_al * self.delta
            for i in range(self.N):
                self.beta += self.eta_be * self.alpha.dot(self.T) ** 2 / 2

            index = self.alpha > 0
            w = (self.alpha * self.T).T.dot(self.X)
            b = (self.T[index] - self.X[index].dot(w)).mean()

            value.append([w, b])

            if (np.power(np.array([k_w - w]), 2) < 0.00001).all() and (k_b - b) ** 2 < 0.00001:
                print(_itr)
                break

            if _itr % 10 == 0 and show:
                self.plot_map(w, b)        
            k_w = w
            k_b = b
        return value
svm = SVM(t=0)
_ = svm.rtin(True)

<はじパタ> おまけ:混合ガウスモデルの検証

前回の記事より。 bython-chogo.hatenablog.com

EMアルゴリズムの実装を行ってみたのですが、アヤメの2次元データではこの混合ガウスモデル、EMアルゴリズムはうまく適用できないね、という結論になりました。2つの山が重なったピーナッツを横にクラスタリングしてほしいのに、縦に割ってしまうようで、初期値をランダムにするとほとんどうまくいかにようでした。

でも、本当でしょうか。私のコードになにか問題があって、うまくクラスタリングができなかった可能性を排除しきれません。いや、むしろこんな超適当コードを信じられないのです。そうだ、私のコードが間違っていてちゃんとEMアルゴリズムの実装ができればうまくいくのではないか。

そんな疑念についてはっきりさせるために積読してあった以下の書籍に書かれているコードを改造して改めて同じデータのプロットをしてみました。

Pythonで学ぶあたらしい統計学の教科書 (AI & TECHNOLOGY)

Pythonで学ぶあたらしい統計学の教科書 (AI & TECHNOLOGY)

その結果が、こ・れ・だ

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

やっぱり縦に割ってるじゃねーか!

というわけでEMアルゴの実力はほかのデータセットで確認をしたほうがよさそうですね。次回ははじパタの11章ですが、もしかしたらサポートベクタマシンのほうに浮気するかもしれません。

<はじパタ> 10.4 確率モデルによるクラスタリング

はい、ここまで10章の教師なし学習のクラスタリングをやってまいりましたが、10章も最後となりました。とりあえず一章をきちんと走り抜けられてよかったです。 で、ここは前回までのK-平均法や融合法と違い、各要素がそれぞれのクラスに入る確率を求める確率モデルによるクラスタリングをおこないます。キーワードとしては、混合正規分布モデルとEMアルゴリズムとなるようです。わかりやすい説明はこちらのリンクをご参照、とりあえず自分なりにまとめてみます。

混合正規分布モデル

クラスタの分布がそれぞれ正規分布にしたがうと仮定をしてその平均と分散が一番Fitする最尤の値を求めていくものになります。また、ある要素に関してそれぞれクラスタの割合を表現する混合比なる重みづけのためのパラメータもでてきてこの値もいじる必要があります。

EMアルゴリズム

所属しているクラスタを示す隠れ変数があるものとして、その変数の確率モデルの最尤推定値を求めるのに優れているのが、EMアルゴリズムとなるようです。E(Expectation)ステップで隠れ変数に対しての期待値を計算、M(Maximization)ステップでその隠れ変数の最大値を計算するようにします。

はじパタに限らず、いくつかの本にこのEMのステップで何が行われているかの数式や説明があるのですが、いかに示すサイトの説明が図として直感的に何をしているのかわかりやすく書かれていると思いました。こういうの貴重。

www.slideshare.net

ピーナッツを縦に割っちゃう問題

それで、コーディングを走らせてみたのですが、どうしても正解率が上がらないんです。実際にプロットをしてみると以下のように右上のピーナッツを縦にわるようなクラスタリングをしてしまうのです。

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

これは初期値にだいぶよるようですね。平均値のmuの値を意図的に変えてみたところ以下のようにピーナッツを横にわるような「理想的な」クラスタに収束することもあるようです。

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

初期値をランダムに置くようにして確認してみましたが、この理想的な形状に収まるのは10回に1~2回程度のようです。

合わないモデル?

初期値問題に頭を抱えて機械学習のお師匠様のところに相談したところ、「そもそも合わないモデルだったのでは」とのことです。考察としてはアヤメデータ以外でもこのEMアルゴリズムを試してみてうまく当てはまるかを確認してみたほうが良いかもしれないです。やるかはわからないですが。

ソースコード

いつも以上にとっちらかったソースになってしまっていると思います。まあ備忘録なので。

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


iris = datasets.load_iris()
features = iris.data
feature_names = iris.feature_names
targets = iris.target

from collections import Counter

class Iris_PModel:
    def __init__(self):
        self.p = []
        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
        
        self.pi = [1/3, 1/3, 1/3]
        self.means_dists()


    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(
                features[p == t, 2],
                features[p == t, 3],
                marker=marker,
                c=c,
            )
            plt.scatter(
                self.mu[t][0], self.mu[t][1], marker='s', c=c
            )
            plt.xlim(1.9, 4.5)
            plt.ylim(0.9, 7.1)
            plt.xlabel(feature_names[2])
            plt.ylabel(feature_names[3])
            plt.autoscale()
            plt.grid()

        plt.show()
    def randomize(self):
        self.pt = np.array([np.random.randint(3) for i in range(sum(self.num))])
        
    def norm(self, mu, sigma):
        n_mu = np.array(mu)
        n_sig = np.array([[sigma[0], sigma[2]],[sigma[2], sigma[1]]])
        n_det = np.linalg.det(n_sig)
        n_inv = np.linalg.inv(n_sig)
        nv = []
        
        for i in range(150):
            n_c = np.array([self.x[i], self.y[i]]) - n_mu
            nv.append((1/(2*np.pi*np.sqrt(n_det))) * np.exp(- n_c.dot(n_inv).dot(n_c).T/2))
        return(nv)
    
    def means_dists(self, init=True):
        if init:
            self.randomize()
            x_m = sum(self.x)/150
            y_m = sum(self.y)/150
            self.nu = [[np.random.normal(x_m, 1), np.random.normal(y_m, 1)], 
                       [np.random.normal(x_m, 1), np.random.normal(y_m, 1)], 
                       [np.random.normal(x_m, 1), np.random.normal(y_m, 1)]]
            print("mu", self.mu)

            self.v  = [[1.0, 1.0, 0.0], [1.0, 1.0, 0.0], [1.0, 1.0, 0.0]]

            
    def m_step(self, gamma):

        n_k = [sum([gamma[i][j] for i in range(150)]) for j in range(3)]

        self.mu = [[sum([self.x[i]*gamma[i][j] for i in range(150)])/n_k[j], 
                    sum([self.y[i]*gamma[i][j] for i in range(150)])/n_k[j]] for j in range(3)] 

        self.v = [[sum([(self.x[i]-self.mu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j],
                  sum([(self.y[i]-self.mu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j],
                  sum([(self.x[i]-self.mu[j][0])*(self.y[i]-self.mu[j][1])*gamma[i][j] for i in range(150)])/n_k[j]] for j in range(3)]                  
                  
        self.pi = [n_k[i]/150 for i in range(3)]
   

    def e_step(self, init=False):
        if init:
            self.install(2, 3)

        nv = [self.norm( self.mu[i], self.v[i]) for i in range(3)]
        snv = [nv[0][i]*self.pi[0] + nv[1][i]*self.pi[1] + nv[2][i]*self.pi[2] for i in range(150)]
       
        prop = [[nv[0][i]*self.pi[0]/snv[i], nv[1][i]*self.pi[1]/snv[i], nv[2][i]*self.pi[2]/snv[i]] for i in range(150)]
        return prop

    def process(self):
        gamma = self.e_step(True)
        self.m_step(gamma)
        
        pre_result = []
        for i in range(31):
            result = []
            for g in gamma:
                result.append(g.index(max(g)))
            self.pt = np.array(result)

            max_int = [0, 0, 0]


            if pre_result == result:
                print(self.pi)
                self.plot_map()
                break

            pre_result = list(result)
            gamma = self.e_step(False)
            self.m_step(gamma)
        

krs = Iris_PModel()
krs.process()

<はじパタ> おまけ:共分散行列のコーディングに苦戦した話

今、はじパタのコーディングをせっせとやってはブログに上げているという作業を続けていますが、いかんせん付け焼刃の知識のために基本的と思われるところで引っかかってしまうことも、多く、今後の備忘録的にも失敗とリカバリの記録を残してみようと思います。

現象:NaNエラーが発生

10.4 「確率モデルによるクラスタリング」をやっていくと、EMモデルの繰り返しによって数値がNaNとなってしまう現象が発生。ひとつNaNがでてくるとその数値を使った値もNaNとなるので真賀田博士もびっくりのすべてがNaNになってしまうのでした。以下がその時のログ。

MU: [[4.2591681450807695, 1.3885253632923946], [3.3363395791482486, 1.0264189358663134], [3.6851362069583895, 1.1830943936011733]]
V:  [[2.8113007747631147, 2.9874538055910858, 3.042873375542265], [0.51266119363098883, 0.56519089148146107, 0.59171106814975671], [1.1252472995735809, 1.2722812742580845, 1.296157580183261]]
N_DET:  -0.860447181357
N_DET:  -0.0603705511147
N_DET:  -0.24839340451
NV: 
nan nan nan
SNV: 
   nan
   nan
   nan
   nan
   nan

原因究明

NaNとなってしまう原因はどこにあるのかを突き止めてみると、ガウス分布の計算内で分散行列の行列式を出すのですが、そこで負の値になってしまっているようでした。その値の平方根を取るものだから当然実数にはならず、NaNとなってしまったのです。なるほど。 そもそものこの共分散行列をかくにんしたところ、

{\displaystyle \Sigma } \Sigma は、半正定値行列
分散共分散行列 - Wikipedia

半正定値行列というのは、行列式が0か正の値になる行列のこと。ということはそもそも負の値が出てくること時点でおかしいということです。 よし、原因はつかめた。あとは訂正をするだけだ。。。ここからが長かった。

行列式のコーディングのチェック

問題となるのは以下の分散を求める部分のコード。意図的にNumpyの行列を使わず、理解のため分散、共分散をそれぞれ計算してみようということで以下のリスト内包表記のコードとなっています。

        self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)]]

最初のsumと次のsumが2次元の要素それぞれの分散値を計算、最後のsumが共分散を計算しています。なんども確認して問題なさそうだったのですが、その結果が以下のように行列式をとると負の値がでてきてしまう。なぜだ。

#分散1、分散2、共分散、行列式(ad-bc)
     2.81130077476 2.98745380559 3.04287337554 -0.860447181357
     0.512661193631 0.565190891481 0.59171106815 -0.0603705511147
     1.12524729957 1.27228127426 1.29615758018 -0.24839340451

入力値とかも確認してみましたがそもそも共分散行列の行列式が負になること自体が異常なので、やはりこのコードがおかしいんですよね。 で、最終的に別の本に書いてあるEMアルゴリズムをコーディングしてみて、その行列要素を出力してみることにしました。すると。。。

V: # Result of ML Book's code 
     2.69263977682 0.640353547689 1.26287084435 0.129398664238
     2.95622008252 0.425917289351 1.08691249439 0.0777264738108
     2.85750222648 0.516449549181 1.19418340666 0.0496817279099

V: # Result of my code
     2.69263977682 2.95622008252 2.85750222648 -0.205283191105
     0.640353547689 0.425917289351 0.516449549181 0.00601751040866
     1.26287084435 1.08691249439 1.19418340666 -0.0534439092193

おおぅ、行列の転置が起きてしまっているようですね。なるほど、リスト内包表記の組み合わせが正しくないようです。基本的な考え方は間違っていないかったですが、こうした順序の認識がまちがっていて、かつ訂正することも難しい(なにが間違っているか理解できない)状況に陥っていました。

訂正

これを踏まえて以下のように訂正をおこなったところ、出力からNaNが出てこなくなりました。

#Corrected
        self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j],
                  sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j],
                  sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j]] for j in range(3)]

#Wrong
        self.v = [[sum([(self.x[i]-self.nu[j][0])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.y[i]-self.nu[j][1])**2*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)],
                  [sum([(self.x[i]-self.nu[j][0])*(self.y[i]-self.nu[j][1])*gamma[i][j] for i in range(150)])/n_k[j] for j in range(3)]]
 

反省

基礎知識がたりないですね。でもまぁ、もともと基礎知識が薄い中でこうしてコケることは織り込み済みで、コケながら覚えていくしかないかと思ってます。今回学んだことをまとめるといかになります。

  • 共分散行列は半正定値行列。行列式が負の値になることはない。
  • リスト内包表記は便利だけど順番を間違えやすい。転置状態に注意。

それでは次回はようやく、「確率モデルによるクラスタリング」にいけるかと思います。

f:id:bython-chogo:20180410163358j:plain

<はじパタ> 10.3 階層型クラスタリング (融合型)

K-平均法(10.2)を扱った前回に引き続きまして、クラスタリングの10章。今回は階層型、融合型というような表現がされている分類方法についてみていきたいと思います。

階層型クラスタリング

この融合型の方法は「定義」により距離的に近いとされる2つのクラスを一つのクラスとして融合することを繰り返していきます。今回は引き続きアヤメのデータを使いましたが、この場合はクラスは3つとなりますので、融合をすすめていき、最終的に3つのクラスに収束したら終わりとしました。

融合するための距離の4つの「定義」

融合にはクラス間の距離的な近さを計測します。はじパタでは距離の定義について以下の4つを紹介しています。

  • 単連結法

2クラスに属する要素通しで最も近いもの通しをそのクラス間の近さと定義して、その値が小さいもの通しを融合していく方法

  • 完全連結法

2クラスに属する要素通しで最も遠いもの通しをそのクラス間の近さと定義して、その値が小さいもの通しを融合していく方法

  • 群平均法

クラスタ内のすべての要素の平均をとり、各クラスの平均値間の距離をクラスの近さと定義する。その値が小さいもの通しを融合していく方法

  • ウォード法 クラスタを融合した場合の平均値が、もとのそれぞれのクラスタ内の平均値から変化する量をクラスタの近さと定義する。その値の小さいもの通しを融合していく方法

コーディングの際の考慮すべき点

アヤメデータは150の要素から成り立つのでそれぞれの要素の距離を示す行列を持つ必要があり、150 x 150のマトリックスとなります。さらにそのうち、上記の定義により近いと判断された2点が同じクラスタとされていきます。要素に番号を振り、融合された要素をリスト化していく形で、マトリックスが縮小(=クラスタ(リスト)内の要素が増えていく)ようにして、最終的に3つのクラスタになった時点で止めるようにしました。 最小値をもとめるコードの大部分は使いまわしが可能で、変更すべき点は、マトリックスのマージのプロセスのみとなりました。

実行結果

アヤメデータ4つの特性をそれぞれ2つずつ、計6通りのセットで確認しています。評価については10.2のK-平均法と同じく正解率とF値を出すようにしています。詳細はこちら

 ===  sepal length (cm)  ===  sepal width (cm)  === 
 ***  Single Linkage Method
   * correctness --  0.35333333333333333
   * F-Value     --  0.502442878626
 ***  Complete Linkage Method
   * correctness --  0.52
   * F-Value     --  0.518563070642
 ***  Group Average Method
   * correctness --  0.7133333333333334
   * F-Value     --  0.741125961406
 ***  Ward Method
   * correctness --  0.78
   * F-Value     --  0.73182865776
 ===  sepal length (cm)  ===  petal length (cm)  === 
 ***  Single Linkage Method
   * correctness --  0.6733333333333333
   * F-Value     --  0.763119756663
 ***  Complete Linkage Method
   * correctness --  0.84
   * F-Value     --  0.825008207203
 ***  Group Average Method
   * correctness --  0.7466666666666667
   * F-Value     --  0.758906000869
 ***  Ward Method
   * correctness --  0.84
   * F-Value     --  0.825008207203
 ===  sepal length (cm)  ===  petal width (cm)  === 
 ***  Single Linkage Method
   * correctness --  0.34
   * F-Value     --  0.501625617698
 ***  Complete Linkage Method
   * correctness --  0.5866666666666667
   * F-Value     --  0.648999234975
 ***  Group Average Method
   * correctness --  0.7066666666666667
   * F-Value     --  0.745888252968
 ***  Ward Method
   * correctness --  0.8266666666666667
   * F-Value     --  0.806976194904
 ===  sepal width (cm)  ===  petal length (cm)  === 
 ***  Single Linkage Method
   * correctness --  0.66
   * F-Value     --  0.760564623487
 ***  Complete Linkage Method
   * correctness --  0.94
   * F-Value     --  0.941265889017
 ***  Group Average Method
   * correctness --  0.8733333333333333
   * F-Value     --  0.863216837848
 ***  Ward Method
   * correctness --  0.8666666666666667
   * F-Value     --  0.854987973338
 ===  sepal width (cm)  ===  petal width (cm)  === 
 ***  Single Linkage Method
   * correctness --  0.66
   * F-Value     --  0.760564623487
 ***  Complete Linkage Method
   * correctness --  0.8066666666666666
   * F-Value     --  0.738999798162
 ***  Group Average Method
   * correctness --  0.66
   * F-Value     --  0.760564623487
 ***  Ward Method
   * correctness --  0.84
   * F-Value     --  0.82020552181
 ===  petal length (cm)  ===  petal width (cm)  === 
 ***  Single Linkage Method
   * correctness --  0.6733333333333333
   * F-Value     --  0.763119756663
 ***  Complete Linkage Method
   * correctness --  0.96
   * F-Value     --  0.960162654599
 ***  Group Average Method
   * correctness --  0.96
   * F-Value     --  0.960162654599
 ***  Ward Method
   * correctness --  0.8933333333333333
   * F-Value     --  0.887000529608

考察

全体的な傾向としては単連結法が比較して正解率が良くなく、ほかの方法はケースバイケースでしょうか。そのなかでもウォード法は安定して高い値がでているように見えます。

K-平均法の時と同じようにすべての特性を扱った4次元データでの精度を確認してみましたが、やはり特徴点が顕著となる2点を使ったデータのほうが精度が高いという結果になりました。

また、K-平均法のほうが精度は高いです。結構苦労してコーディングしたのに悲しい限り。アヤメデータ以外でのK-平均法よりよくなるケースがあるのか気になるところではあります。

ソースコード

いかJupyter notebook上で動くソースです。

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

class ClastaringFusion():
    def __init__(self):
        self.x = []
        self.y = []
        self.p = []
        self.t = []
        self.pt = []
        self.num = [50, 50, 50]
        iris = datasets.load_iris()
        
        self.feature_names = iris.feature_names
        features = iris.data
        targets = iris.target
        self.dmatrix = np.array([])

    def install_val(self,a=1,b=2):
        self.p = [a, b]
        self.x = features.T[self.p[0]]
        self.y = features.T[self.p[1]]
        self.t = targets
        self.distance()

        
    def distance(self, x=None, y=None):
        if x==None: x = self.x

        self.dmatrix = np.array([[np.sqrt((self.x[i] - self.x[j])**2 + (self.y[i] - self.y[j])**2 ) for i in range(150)] for j in range(150)]) 
        
        #return np.array([[np.sqrt((i - j)**2 + (i - j)**2) for i in x] for j in y])          
    def comp_fusion(self, list1, list2, fm=0):
        def ave(lst):
            x_lst = [self.x[i] for i in lst]
            y_lst = [self.y[i] for i in lst]

            m_x = sum(x_lst)/len(x_lst)
            m_y = sum(y_lst)/len(y_lst)


            return sum([(x_lst[i] - m_x)**2 + (y_lst[i] - m_y)**2 for i in range(len(lst))])
            
        if fm ==0:
            v = min([min([self.dmatrix[l1, l2] for l2 in list2]) for l1 in list1])
        elif fm==1:
            v =  max([max([self.dmatrix[l1, l2] for l2 in list2]) for l1 in list1])
        elif fm==2:
            v = sum([sum([self.dmatrix[l1, l2] for l2 in list2]) for l1 in list1]) /(len(list1) * len(list2))
        elif fm==3:
            v = ave(list1+list2) - (ave(list1) + ave(list2))    
        else:
            print("Error: Invalid number for Fusion Method")
            sys.exit(1)
            
        return v
    
    def fusion(self, fm=0, tree=None):
        '''
        fm (fusion_method)
        0: Single Linkage Method
        1: Complete Linkage Method
        2: Group Average Method
        3: Ward Method
        
        '''
        if tree==None:
            tree = [[i] for i in range(150)]
        dmtx = np.array([[self.comp_fusion(i, j, fm) for i in tree] for j in tree])

        for i in range(dmtx.shape[0]):
            dmtx[i, i] = 100
        
        val = np.argmin(dmtx)
        #print(val)
        #print(dmtx.shape)
        #print("D", dmtx[val//dmtx.shape[0], val%dmtx.shape[0]])
        minval = np.array(np.where(dmtx==dmtx[val//dmtx.shape[0], val%dmtx.shape[0]]))

        ctree = []
        checked = []
        lgth = minval.shape[1]
        #print(minval)
        #print("lgth: ", lgth)
        for i in range(lgth): 
            flag = True
            one = tree[minval[0][i]]
            two = tree[minval[1][i]]
            checked += (one + two)
            inc_num = []
            
            #print("one+two", one, two)
            #print(" ctree ", ctree)
            if minval[0][i] != minval[1][i]:

                for j, t in enumerate(ctree):
                    if set(t) & set(one) or set(t) & set(two):
                        inc_num.append(j)
                #print(" ", inc_num, " ==== ", one, two)
                #print(" ", ctree)
                if len(inc_num) == 0:
                    ctree.append(one+two)
                else:
                    ctree[inc_num[0]] = list(set(one+two+ctree[inc_num[0]]))
                    if len(inc_num) > 1:
                        for n in inc_num:
                            ctree[inc_num[0]] = list(set(ctree[inc_num[0]]+ctree[n]))
                            
                        for k, n in enumerate(inc_num[::-1]):
                            if k == len(inc_num)-1:
                                break
                            ctree.pop(n)
                    
          
        for i in tree:
            if not set(i) & set(checked):
                ctree.append(i)
        #print(ctree)
        test = []

        for t in ctree:
            test += t
        
        missing =  (set(range(150)) - set(test))
        tomatch =  (set(test) - set(range(150)))
        #print("Missing", missing)
        #print("Too Match", tomatch)
        return ctree  
    def show_all_result(self):
        fm_string = [
            "Single Linkage Method",
            "Complete Linkage Method",
            "Group Average Method",
            "Ward Method"
        ]
        for i, (a, b) in enumerate(itertools.combinations(range(4), 2)):
            
            print(" === ", self.feature_names[a], " === ", self.feature_names[b], " === ")
            self.install_val(a, b)
            #self.plot_map()

            for fm in range(4):
                tree = cf.fusion(fm)
                c = 0
                while len(tree) > 3:
                    #print("======")
                    tree = cf.fusion(fm, tree)
                    #print(len(tree))
                    #if len(tree) < 4:
                    #    break

                    c += 1
                crct, fv = cf.correct_rate(tree, False)
                print(" *** ", fm_string[fm])
                print("   * correctness -- ", crct)
                print("   * F-Value     -- ", fv)
             
        
    def correct_rate(self, tree, verbose=True):

        trees = []
        crcts = 0.0
        fvs = 0.0
        

        for p in list(itertools.permutations(range(3))):

            pt_tree = []
            for i in range(150):
                for j in range(3):
                    if i in tree[j]:
                        pt_tree.append(p.index(j))
                        break
                        
            trees.append(pt_tree)
            #print("tree", pt_tree)
                
            cnt = np.array([[Counter(pt_tree[i*50:i*50+50])[j] for i in range(3)] for j in range(3)])
            odr = [np.argmax([Counter(pt_tree[j*50:j*50+50])[i] for i in range(3)]) for j in range(3)]
            crct = []
            crctnum = 0

            for i in range(150):
                if i // 50 == pt_tree[i]:
                    crctnum += 1

            fv = 1.0

            for i in range(3):

                x = np.count_nonzero(pt_tree[i*50:i*50+50] == odr[i])
                z = np.count_nonzero(pt_tree[(i+1)%3*50:(i+1)%3*50+50] == odr[i]) + \
                    np.count_nonzero(pt_tree[(i+2)%3*50:(i+2)%3*50+50] == odr[i])
                rec = x/50
                pre = 50/(50+z)
                fv *= 2*rec*pre/(rec+pre)

            if crctnum > crcts:
                fvs = fv
                crcts = crctnum
            #print("crct", crctnum)
        if verbose:
            print("CORRECTNESS: ", crcts/150)
            print("F-Value    : ", np.cbrt(fvs))
        return crcts/150, np.cbrt(fvs)
        
cf = ClastaringFusion()
cf.install_val(1, 2)
cf.distance()
cf.show_all_result()

関連記事

bython-chogo.hatenablog.com

bython-chogo.hatenablog.com

<はじパタ> 10.2 K-平均法 おまけ

今回は以下の記事のおまけになります。

bython-chogo.hatenablog.com

K-平均法はランダムにK個クラスタした平均値をとって、その平均値に近い距離にいる要素をそのクラスタに分類する、という試行を収束するまで繰り返すものでした。

じゃあ実際にその平均値は各試行ごとにどのように変化するのかというのをアニメーションでプロットしてみました。 それが以下です。四角いポイントが各色の平均値となり、繰り返しによるクラスタの変化を各ステップごとにアニメーションで示しています。

f:id:bython-chogo:20180403235608g:plain

なるほど、最初にランダムにアサインされたデータから各テリトリーごとに収束していく様子が見やすくなりました。これは面白いですね。

ただ、ランダムでアクセスするので、以下のようにちょっとおかしな収束してしまうこともあるようです。特殊なパターンの収束をさけるために、試行は一回ではなく複数回やるのがよさそうです。

f:id:bython-chogo:20180403235600g:plain

さて、次回は10.3融合法です。この方法もアヤメデータを使った分類をしてみて、その正確さを評価してみたいと思います。