「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

統計的機械学習の数理100問(20/20)

 いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

を用いることにする。

9. 教師なし学習

 教師データが無い中で説明変数\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n\in\mathbb{R}^pのみからそれら標本間の関係、p変数間の関係を学習する問題を扱う。

9.1 K-meansクラスタリング

 クラスタリングp変数を持つn個の標本\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n\in\mathbb{R}^pを交わりの無いK個の集合(クラスター)に分割する処理を指す。このうちK-meansクラスタリングは事前にK個の値を定め、最初に1,\cdots,Kのいずれかをnこの標本それぞれにランダムに割り当てて、以下のステップを繰り返す:

   (1) クラスタk=1,\cdots,Kのそれぞれに含まれるサンプルの平均ベクトルを計算する。
   (2) n個の標本それぞれにKクラスターの中で平均ベクトルが最も近いクラスターを当てはめる。

ここでクラスターはp次元ベクトルの集合で、その平均はそれらの中心であることを意味する。また2段階目における近さは\mathcal{L}^2-ノルム


\begin{aligned}
\|\boldsymbol{a}-\boldsymbol{b}\|=\displaystyle{\sqrt{\sum_{i=1}^{p}(a_i-b_i)^2}}
\end{aligned}

で測る。n\lt kである場合、所属するサンプルが存在しないクラスターが生じる場合がある。このときはそのクラスターを無視する。
 K-meansクラスタリングを実行している間、毎回の更新でスコア 


\begin{aligned}
S=\displaystyle{\sum_{k=1}^{K}\min_{\boldsymbol{z}_k\in\mathbb{R}^p}\left(\sum_{i\in C_k}\|\boldsymbol{x}_i-\boldsymbol{z}_k\|^2\right)}
\end{aligned}


は増加しない。またK-meansクラスタリングの結果はランダムに選ばれた初期状態のクラスタに依存する。このことはK-meansクラスタリングを適用しても最適解が得られる保証がなく、初期値を変更して何度も実行し、その中でスコアの良いクラスタリングを選ぶ必要があることを意味する。


 \boldsymbol{x}_1={}^{t}(x_{1,1},\cdots,x_{1,n}),\cdots,\boldsymbol{x}_n={}^{t}(x_{n,1},\cdots,x_{n,p})および各クラスタk=1,\cdots,Kが含むサンプルの添字集合をC_k、各クラスタの中心を\bar{\boldsymbol{x}}_1={}^{t}(\bar{x}_{1,1},\cdots,\bar{x}_{1,p}),\cdots,\bar{\boldsymbol{x}}_n={}^{t}(\bar{x}_{n,1},\cdots,\bar{x}_{n,p})とおけば、



\begin{aligned}
\displaystyle{\frac{1}{|C_k|}\sum_{i\in C_k}\sum_{i^{\prime}\in C_k}\sum_{j=1}^{p}(x_{i,j}-x_{i^{\prime},j})^2}=2\displaystyle{\sum_{i\in C_k}\sum_{j=1}^{p}(x_{i,j}-x_{i^{\prime},j})^2}
\end{aligned}


が成立する。このため、K-meansクラスタリングクラスター内のサンプルに関するすべての対の距離の二乗和を最小にすることを目的としている。

9.2 階層的クラスタリング

 K-meansクラスタリングに並んで頻用されるのが階層的クラスタリングである。
 最初に1サンプルのみを含むクラスターをn個用意し、ある基準に基づき近いクラスターを結合していき、クラスター数を2まで減らす。事前にクラスター数Kを決めずにすべての1\leq K\leq nでのクラスタリングを同時に見出すことができる。サンプル間の距離は普通、\mathcal{L}^2-ノルムを用いる。しかし複数サンプルを含んでいるクラスター同士の距離を定義する必要がある。

9.3 主成分分析

 主成分分析はp\leq nとして、ある行列X\in\mathbb{R}^{n\times p}が与えられたときに、\left\|X\phi_1\right\|を最大にするような\phi_1\in\mathbb{R}^n,\phi_iと直交して\left\|X\phi_2\right\|を最大にする\phi_2\in\mathbb{R}^n,\cdots,\left\|X\phi_p\right\|を最大にする\phi_p\in\mathbb{R}^nを求める操作である。ただし\|\phi_1\|=\cdots=\|\phi_\|=1を仮定する。
 主成分分析の主要な目的は、\phi_1から\phi_m,1\leq m\leq pまでの成分を用いて行列Xに関する情報を要約することである。このとき



\begin{aligned}
{}^{t}XX\phi_1=\mu_1\phi_i
\end{aligned}


を満たすような定数\mu_1が存在する。そのような\mu_1は複数あり得るが、


\begin{aligned}
\|X\phi_1\|^2={}^{t}\phi_1{}^{t}XX\phi_1=\mu_1\|\phi_1\|^2=\mu_1
\end{aligned}

を最大にするため、最大の固有値を選ぶ必要がある。
 また\phi_2もある\mu_2について


\begin{aligned}
{}^{t}XX\phi_2=\mu_2\phi_2
\end{aligned}

を満たさなければならない。したがって{}^{t}XX固有ベクトルになる。しかし\mu_1\geq\mu_2であって\phi_1と直交する必要がある。すなわち\mu_1=\mu_2\phi_1,\phi_2が同じ固有空間に属するか、\mu_2\mu_1以外で最大の固有値である必要がある。またそれらは非負定値行列の固有値であるため、非負である。
 実際の主成分分析の定式化では、\displaystyle{\frac{1}{n}{}^{t}XX}とおいて、\displaystyle{\frac{\mu_1}{n}},\cdots,\displaystyle{\frac{\mu_m}{n}}を改めて\lambda_1,\cdots,\lambda_mとおき、



\begin{aligned}
{}^{t}XX\phi_1=\mu_1\phi_i
\end{aligned}



\begin{aligned}
\Sigma\phi_1=\lambda_1\phi_1
\end{aligned}


と書くことが多い。ここで\Sigmaは標本分散共分散行列であり、\lambda_1\geq\cdots\geq\lambda_m\geq0はそれぞれ第1主成分\phi_1,\cdots,m主成分\phi_mの分散と呼ばれる。


 主成分分析には別の解釈がある。行列X\in\mathbb{R}^{n\times p}の第i行、\Phi\in\mathbb{R}^{p\times m}を各列[tex:\phi_1,\cdots,\phi_mの大きさが1で直交するベクトルだとする。このとき\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n\phi_1,\cdots,\phi_mへの射影\boldsymbol{z}_1=\boldsymbol{x}_1\Phi,\cdots,\boldsymbol{z}_n=\boldsymbol{x}_n\Phi\in\mathbb{R}^mが得られる。
 \boldsymbol{z}_1,\cdots,\boldsymbol{z}_nの右から{}^{t}\Phiを掛けて\boldsymbol{x}_1,\cdots,\boldsymbol{x}_nに近いベクトルが復元できたかどうかを



\begin{aligned}
L=\displaystyle{\sum_{i=1}^{n}\left\|\boldsymbol{x}_i-\boldsymbol{x}_i\Phi{}^{t}\Phi\right\|^2}
\end{aligned}


で評価する。もしm=pであるならば、この値は0になる。主成分分析はこの値を最小にするような\phi_1,\cdots,\phi_mを求める問題と見ることができる。

9.4 Rによるシミュレーション

9.4.1 k-meansクラスタリング
############################
### K-meanクラスタリング ###
############################

library("ggplot2")

k_means <- function(X, K, iteration = 20) {
  n <- nrow(X)
  p <- ncol(X)
  center <- array(dim = c(K, p))
  y <- sample(1:K, n, replace = T)
  for (h in 1:iteration) {
    for (k in 1:K) {
      if (sum(y[] == k) == 0) {
        # sum(y[] == k)でy[i] = kなるiの個数を意味する。
        # サンプルを含まないクラスタは消滅
        center[k, ] <- Inf
      } else {
        for (j in 1:p)
          center[k, j] <- mean(X[y[] == k, j])
      }
    }
    for (i in 1:n) {
      S_min <- Inf
      for (k in 1:K) {
        S <- sum((X[i, ] - center[k, ]) ^ 2)
        if (S < S_min) {
          S_min <- S
          y[i] <- k
        }
      }
    }
  }
  return(list(clusters = y))
}

# データ生成
p <- 2
n <- 1000
X <- matrix(rnorm(p * n), nrow = n, ncol = p)

# 各サンプルのクラスタを得る
y <- k_means(X, 5)$clusters

df_out <- data.frame(x = X[,1],
                     y = X[,2],
                     lab = y)
df_out$lab <- factor(df_out$lab,level = 1:5)

# クラスタごとに色を変えて,点を描く
g <- ggplot(df_out,mapping = aes(x = x, y = y, color = lab))
g <- g + geom_point()
g <- g + labs(x = "第1成分", y = "第2成分")
g <- g + theme_classic()
g <- g + theme(plot.title = element_text(hjust = 0.5),axis.text = element_text(size = 10),
               axis.title.x = element_text(size = 12),
               axis.title.y = element_text(size = 12),
               legend.position = "none",
               legend.title = element_text(size = 10),
               legend.text = element_text(size = 10))
g <- g + geom_vline(xintercept = 0, linetype = 3)
g <- g + geom_hline(yintercept = 0, linetype = 3)

plot(g)

9.5 Pythonによるシミュレーション

9.5.1 k-meansクラスタリング
##############
### kMeans ###
##############

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn #正規乱数


def k_means(X,K,iteration=20):
    n,p=X.shape
    center=np.zeros((K,p))
    y=np.random.choice(K,n,replace=True)
    scores=[]
    for h in range(iteration):
        for k in range(K):
            if np.sum(y==k)==0:
                center[k,0]=np.inf
            else:
                for j in range(p):
                    center[k,j]=np.mean(X[y==k,j])
        S_total=0
        for i in range(n):
            S_min=np.inf
            for k in range(K):
                S=np.sum((X[i,]-center[k,])**2)
                if S<S_min:
                    S_min=S
                    y[i]=k
            S_total+=S_min
        scores.append(S_total)
    return {'clusters':y, 'scores':scores}

# テスト
n=500
K=5
p=2
X=randn(n,p)
y=k_means(X,5)['clusters']

#クラスタごとに色を変えて,点を書く
plt.scatter(X[:,0],X[:,1],c=y)
plt.xlabel("第1成分")
plt.ylabel("第2成分")
9.5.2 階層的クラスタリング
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import japanize_matplotlib
import scipy
from scipy import stats
from numpy.random import randn #正規乱数

import copy

def dist_complete(x,y):
    r=x.shape[0]
    s=y.shape[0]
    dist_max=0
    for i in range(r):
        for j in range(s):
            d=np.linalg.norm(x[i,]-y[j,])
            if d>dist_max:
                dist_max=d
    return dist_max

def dist_single(x,y):
    r=x.shape[0]
    s=y.shape[0]
    dist_min=np.inf
    for i in range(r):
        for j in range(s):
            d=np.linalg.norm(x[i,]-y[j,])
            if d<dist_min:
                dist_min=d
    return dist_min

def dist_centroid(x,y):
    r=x.shape[0]
    s=y.shape[0]
    x_bar=0
    for i in range(r):
        x_bar=x_bar+x[i,]
    x_bar=x_bar/r
    y_bar=0
    for i in range(s):
        y_bar=y_bar+y[i,]
    y_bar=y_bar/s
    return (np.linalg.norm(x_bar-y_bar))

def dist_average(x,y):
    r=x.shape[0]
    s=y.shape[0]
    S=0
    for i in range(r):
        for j in range(s):
            S=S+np.linalg.norm(x[i,]-y[j,])
    return(S/r/s)

def hc(X,dd="complete"):
    n=X.shape[0]
    index=[[i] for i in range(n)]
    cluster=[[] for i in range(n-1)]
    for k in range(n,1,-1):
        #index_2=[]
        dist_min=np.inf
        for i in range(k-1):
            for j in range(i+1,k):
                i_0=index[i];j_0=index[j]
                if dd == "complete":
                    d=dist_complete(X[i_0,],X[j_0,])
                elif dd == "single":
                    d=dist_single(X[i_0,],X[j_0,])
                elif dd == "centroid":
                    d=dist_centroid(X[i_0,],X[j_0,])
                elif dd == "average":
                    d=dist_average(X[i_0,],X[j_0,])
                if d<dist_min:
                    dist_min=d
                    i_1=i    #結合される側のlistのindex
                    j_1=j    #新たに結合するlistのindex
        index[i_1].extend(index[j_1])  #追加する
        if j_1 < k:               #追加したindexの後ろを1つ前に詰める
            for h in range(j_1+1,k,1):
                index[h-1]=index[h]
        index2=copy.deepcopy(index[0:(k-1)])  #indexのまま使うと, 毎回書き換わってしまうため
        cluster[k-2].extend(index2)
    return cluster  #下から結果を見ると, 1つずつ結合が起こっていることがわかる
n=200;p=2
X=randn(n,p)
cluster=hc(X,"complete")
K=[2,4,6,8] #クラスタ数が3,5,7,9

for i in range(4):
    grp=cluster[K[i]]  #全体の結果から, クラスタ数がK[i]のときの結果を取り出す
    plt.subplot(2,2,i+1)
    for k in range(len(grp)):
        x=X[grp[k],0]
        y=X[grp[k],1]
        plt.scatter(x,y,s=5)
    plt.text(2,2,"K={}".format(K[i]+1),fontsize=12)

n=100;p=2;K=7
X=randn(n,p)
i=1
for d in ["complete","single","centroid","average"]:
    cluster=hc(X,dd=d)
    plt.subplot(2,2,i)
    i=i+1
    grp=cluster[K-1]
    for k in range(K):
        x=X[grp[k],0]
        y=X[grp[k],1]
        plt.scatter(x,y,s=5)
    plt.text(-2,2.1,"{}".format(d),fontsize=12)
9.5.3 主成分分析
def pca(X):
    n,p=X.shape
    center=np.average(X,0)
    X=X-center #列ごとに中心化
    Sigma=X.T@X/n
    lam,phi= np.linalg.eig(Sigma)  #固有値, 固有ベクトル
    index = np.argsort(-lam)  #降順にソート
    lam = lam[index]
    phi = phi[:,index]
    return {'lam':lam,'vectors':phi,'centers':center}

X=randn(100,5)
res=pca(X)
res['lam']  

res['lam']/np.sum(res['lam'])  #各主成分の寄与率, 下のevrと同じ
res['vectors']
res['centers']

# 
from sklearn.decomposition import PCA
pca=PCA()
pca.fit(X) #実行
score=pca.fit_transform(X) #主成分得点(行:n 列:主成分)
score[0:5,]
pca.components_

pca.mean_ 

evr=pca.explained_variance_ratio_  #各主成分が全体のどれだけ説明しているか
evr

plt.plot(np.arange(1,6),evr)
plt.scatter(np.arange(1,6),evr)
plt.xticks(np.arange(1,6))
plt.ylim(0,1)
plt.xlabel("主成分")
plt.ylabel("寄与率")

##################

n=100;a=0.7
b=np.sqrt(1-a**2)
u=randn(n);v=randn(n)
x=u
y=u*a+v*b
plt.scatter(x,y)
plt.xlim(-4,4)
plt.ylim(-4,4)

D=np.concatenate((x.reshape(-1,1),y.reshape(-1,1)),1)
pca.fit(D)

T=pca.components_   
T[0,1]/T[0,0]*T[1,1]/T[1,0]   #主成分ベクトル(主成分負荷量)が直交している

def f_1(x):
    y=T[0,1]/T[0,0]*x
    return y
def f_2(x):
    y=T[1,1]/T[1,0]*x
    return y
x_seq=np.arange(-4,4,0.5)
plt.scatter(x,y,c="black")
plt.xlim(-4,4)
plt.ylim(-4,4)
plt.plot(x_seq,f_1(x_seq))
plt.plot(x_seq,f_2(x_seq))
plt.gca().set_aspect('equal', adjustable='box')
プライバシーポリシー お問い合わせ