＜はじパタ＞ おまけ：混合ガウスモデルの検証

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

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

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

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

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

その結果が、こ・れ・だ

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

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

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

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

混合正規分布モデル

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

EMアルゴリズム

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

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

数式を使わずイメージで理解するEMアルゴリズム from 裕樹 奥田

www.slideshare.net

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

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

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

初期値をランダムに置くようにして確認してみましたが、この理想的な形状に収まるのは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

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

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

問題となるのは以下の分散を求める部分のコード。意図的に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)]]
 

反省

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

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

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

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

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

階層型クラスタリング

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

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

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

• 単連結法

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

• 完全連結法

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

• 群平均法

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

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

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

アヤメデータは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個クラスタした平均値をとって、その平均値に近い距離にいる要素をそのクラスタに分類する、という試行を収束するまで繰り返すものでした。

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

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

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

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

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

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

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

アヤメデータ

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

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

つまり、今回の問いはこの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-平均法の結果

• 実際のデータ

• 評価値

 ===  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値はピークを越えて少しだけ下がるところで収束してしまうようです。実際は回答とすべき教師データはないので、収束まつしか判断のしようがなさそうです。

４つの特性から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といったあたりに収まるでしょうか。

４つの特性を使った評価

では、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) 

CTS LabsのAMDプロセッサに対する脆弱性のPaperをまとめてみた

3月13日に公表されたAMDの脆弱性について、イスラエルのセキュリティ会社、CTS Labsが出しているWhite Paperをまとめてみた。

AMDFLAWS

POC（Proof of Concept: コンセプトの 実効性検証）の開示もまだのようで具体的な手法の情報も乏しく、現時点で実現可能性に関しては疑問がのこるものになる。

いくつかのニュースサイトで記事になっているように、現時点での情報ではこの攻撃による深刻性については不明確だ。すでに強い権限を持った状態が前提の作業で、前段階のセキュリティが破られいることは条件になりそうだ。また、この手の脆弱性にはCert/CCから発行されるCVEナンバーがついて、公表されるのが定例だが、この会社はそのような公開プロセスを飛ばして公表に踏み切っているというため、お行儀がよくないと批判も受けているとのこと。 なんにせよ、最終的な結論は詳細情報のアップデートが必要そうだ。

CTSの以下の記述のとおり、現在、対応パッチや緩和策のためにAMD/Microsoftなどの関係企業でシェアしているとのこと。（3月15日現在）

Q. Doesn't this publication put users at risk?
A. No. All technical details that could be used to reproduce the vulnerabilities have been redacted from this publication. CTS has shared this information with AMD, Microsoft, and a small number of companies that could produce patches and mitigations.

以下、公式サイトおよび論文のまとめとなる。技術知識が深くない私が、Paperを読んで理解できた上っ面のまとめとなること、勘違い、知識不足をご容赦いただきたい。

脆弱性概要

CTS Labsによると今回公表されている脆弱性は大きく4つのクラス(MASTERKEY/RYZENFALL/FALLOUT/CHIMERA)に分けられ、13個の脆弱性があるという。

公式サイトに記載されている各クラスの概況を大雑把にまとめる。

• MASTERKEY: Unauthorized Code Execution and Malware Persistency on AMD Secure Processor

-- AMD Secure Processor firmwareに対する攻撃で、Firmware ベースのセキュリティ機能 (Secure Encrypted Virtualization, Firmware Trusted Platform Module)に対する書き換え、FirmwareベースでのRansomewareの書き込みが可能だという。

• RYZENFALL: Vulnerabilities in Ryzen Secure Processor

-- SMRAM(System Momory RAM), やWindowsが持っているセキュアなメモリ空間にアクセスすることが可能になるため、AMD Secure Processorをコントロールすることができる。MASTERKEY攻撃に必要な不正なBIOSコードの書き込みも行える。

• FALLOUT: Vulnerabilities in EPYC Server Secure Processor

-- SMRAM(System Momory RAM), やWindowsが持っているセキュアなメモリ空間にアクセスすることが可能になるため、AMD Secure Processorをコントロールすることができる。不正なBIOSコードの書き込みも行える。RYZENFALLはSecure OS機構に対する攻撃でClient/Workstation製品に対する攻撃だが、FALLOUTはBoot Loaderチップに対する攻撃で対象はEPYC Server製品となる

• CHIMERA: Backdoors Inside Ryzen Chipset

-- Ryzenなどのプラットフォームで用いられている3rd Partyの隠れたBackdoorを突くもの。DMA(Direct Memory Access)エンジンなどに攻撃をすることで、プラットフォームのIOに対する不正攻撃が可能となる。Key Logger/Network Access/PypassMemory Protectionなどが挙げられている。

Paper まとめ

MASTERKEY: Unauthorized Code Execution and Malware Persistency on AMD Secure Processor

AMDのSecure Bootの機能にMalwareなど不正なコードをしこみ、最終的にSecure ProcessorのKernel Modeの挙動もコントロールが可能とされている。Malwareを仕込めば、Secure Bootを不正にBypassしてSecure Processorの機能を無効化もできる。 通常はSignatureでサインされているものでないとBIOSのアップデートはできないのだが、RYZENFALLもしくはFALLOUTでSystem Management Mode(SMM)という上級権限に入ることにより、それを可能としている。

RYZENFALL: Vulnerabilities in Ryzen Secure Processor

RyzenなどClient/Workstation向けのシステムに導入されているSecure OSはfTPM (Firmware Trusted Platform Moduel)に含まれる認証情報からのスタートとなる。このSecure BootのProcessで利用されるメモリ空間は特に権限が強い Fenced DRAMというところでSecure Processorでしか利用ができない空間になる。 ここに対するRYZENFALLの攻撃はadministrator 権限が必要となり、Secure Processorへのアクセスはベンダから公式に配布されるドライバを使用する。RYZENFALLに成功するとSecure Processorでの実行が認証なしで行えるようになり、ハードウェアで保護されているMemory Regionへのアクセスを可能とする。具体的には以下の攻撃チャネルがある。

• Windows Isolated User Mode and Isolated Kernel Mode (VTL1)
• Secure Management RAM (SMRAM)
• AMD Secure Processor Fenced DRAM

これらを経由することで、 以下のことが実行可能になる。 - Microsoft Virtualization-based Securityをバイパスして Networkの権限を奪うことが可能。 - SMMにMalwareをいれることができOSやHypervisorでの攻撃が可能となる - BIOSのReflashに対するプロテクションの無効化 - VTL1にMalwareを入れることが可能。OS上でのエンドポイントでのSecurity機能の無効化が可能。 - AMD Secure Processor自身にMalwareの注入が可能 - fTPMなどのSecurity機能をバイパスしてのブートアップが可能

FALLOUT: Vulnerabilities in EPYC Server Secure Processor

FALLOUTはEPYCなどのサーバを対象とした攻撃方法。EPYCのSecurity Processorで利用するBoot loaderのコンポーネント内の脆弱性をつくと、Hardware Validated BootやSecure Encrypted VirtualizationなどのSecurityに影響を与えることができる。 RYZENFALLと同じく、administrator権限での実行が必要となり、公式のドライバを使用する。これを利用してHardwareで保護されているMemory Regionへのアクセスすることができる。具体的には以下の攻撃チャネルがある。

• Windows Isolated User Mode and Isolated Kernel Mode (VTL1)
• Secure Management RAM (SMRAM)

これらを経由することで、 以下のことが実行可能になる。 - Microsoft Virtualization-based Securityをバイパスして Networkの権限を奪うことが可能。 - SMMにMalwareをいれることができOSやHypervisorでの攻撃が可能となる - BIOSのReflashに対するプロテクションの無効化 - VTL1にMalwareを入れることが可能。OS上でのエンドポイントでのSecurity機能の無効化が可能。

CHIMERA: Backdoors Inside Ryzen Chipset

AMDのPromontory Chipset内にある、製造工程上などで利用されたとみられる隠れたバックドアを突く方法。2つのバックドアが存在している。一つはChipsetを実行するFirmwareにあり、ASICハードウェアハードの中にある。特に後者はChipset内にあるものなので修正が難しく、ワークアラウンドかRecallしかない。 Chipset内にあるので、提供されているSecurityの機能にも影響を及ぼすことが可能になり、さらにChipsetはUSB, SATA, PCI-EといったI/O周りも制御ができるので、LAN, Wifi, Bluetoothにも影響を及ぼすことができる。彼らの研究ではDMA（Direct Memory Access）に対する攻撃で、OSに影響を与えることができたという。このことにより、以下のような追加の攻撃ができるようになると考えられている。

• Key Logger
• Network Access
• Bypass Memory Protection