統計的機械学習の数理100問（20/20） - 「大人の教養・知識・気付き」を伸ばすブログ

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

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

作者:讓, 鈴木
共立出版

Amazon

を用いることにする。

前回

https://power-of-awareness.com/entry/2022/07/05/120000power-of-awareness.com

前回
9.　教師なし学習

9.　教師なし学習

教師となるデータが無く、 $\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n$ のみで $n$ 個のサンプルの間にある関係や $p$ 変数間の関係を学習する方法を教師無し学習という。

ここではクラスタリングと主成分分析を取り上げる。

クラスタリングは $\boldsymbol{x}_1,\cdots,\boldsymbol{x}_n$ をいくつかのグループに分ける。クラスタ数を事前に与える $K$ - $\mathrm{means}$ クラスタリングとそのような情報を与えない階層的クラスタリングの2つがある。

主成分分析は機械学習に限らず多変量解析などで頻用される。

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

　クラスタリングは $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  #下から結果を見ると, １つずつ結合が起こっていることがわかる
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')

前回

9. 教師なし学習

9.1 K-meansクラスタリング

9.2 階層的クラスタリング

9.3 主成分分析

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

9.4.1 k-meansクラスタリング

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

9.5.1 k-meansクラスタリング

9.5.2 階層的クラスタリング

9.5.3 主成分分析