機械学習レポート

※日本ディープラーニング協会のDeep Learning資格試験（E資格）受験に必要な、協会認定の講座プログラム( https://ai999.careers/rabbit/ )がある。
※本プログラム受講に際し、レポート作成はWebやBlogなどで作成する必要があるとのことで、Qiitaで作成することとした。
※個人的な備忘録兼レポート作成用であり、網羅的にカバーするものではありません。ご了承ください。

機械学習については、線形回帰モデル、非線形回帰モデル、ロジスティック回帰モデル、主成分分析、アルゴリズム、サポートベクターマシーンの各項目についてそれぞれ100文字以上で要点をまとめるのに加え、実装演習結果キャプチャーまたはサマリーと考察が必要とのことである。

以下に各項目ごとに記載する。

1.線形回帰モデル

線形回帰モデルは最適化問題であらゆる分野で広く使用されてきた手法であり、機械学習に限る問題ではないが、機械学習の1手法という整理の仕方がされている。

ある入力から出力を予測する問題を回帰問題といい、入力と出力の関係を1次式で表せるものを線形回帰モデルで扱う。まず、m個の入力値と出力yとの線形関係式を構築する。入力は説明変数とも呼ばれ入力の数がm個の場合は

\boldsymbol{x}=(x_1,x_2,\ldots,x_m)^{T} \in\mathbb{R}^{m}

ただし、$\in\mathbb{R}^{m}$はm次元実数空間、実数m個の入力値を示す。目的変数$y\in\mathbb{R}^{1}$は、実数1次元の出力を意味する。それぞれの入力$x_i$に対して係数（あるいはパラメータ、重みなどと呼ばれる）$w_i$がかかるので、$\boldsymbol{w}$は$\boldsymbol{x}$と同じ次元をもつ。

\boldsymbol{w}=(w_1,w_2,\ldots,w_m)^{T} \in\mathbb{R}^{m}

1つのデータセットにおける入力、出力、係数との関係は以下で表される。

\hat{y}=\boldsymbol{w}^{T}\boldsymbol{x}+w_0

ただし、$\hat{y}$は入力値とともに与えられた既知の出力データ（教師データ）、$w_0$は入力にかからないシフト量とする。

入力と出力のデータ（教師データ）$n$個ある場合は以下のような行列表現ができる。

\boldsymbol{y}=\boldsymbol{Xw}+\boldsymbol{\epsilon}

ただし、$\boldsymbol{y}$はn x 1の出力値を表すベクトル（データ数n個に対応）、$\boldsymbol{X}$は入力値を表すn x mの行列（それぞれの教師データ毎にm個の入力値があり、データ数がn個ある）、$\boldsymbol{w}$は係数を表す(m+1) x 1のベクトル（各入力値に対する重みの数とバルクシフト$w_0$）である。

既知である入力・出力データ（教師データと呼ばれる）から、最適な重み$\boldsymbol{w}$を最小二乗法で求める。

\boldsymbol{\hat{w}}=\arg \min MSE

MSEはMean Square Error（最小二乗誤差の平均値）

MSE=\frac{1}{n}\sum_{i=1}^n{(\hat{y}_i (\boldsymbol{w})-y_i)^2}

また、argminは、最小となるMSEの値そのものではなく、その時の変数$\boldsymbol{w}$を返すという意味である。MSEが最小となるのはMSEを$\boldsymbol{w}$で微分して0になる点であるため、

{\frac{\partial}{\partial\boldsymbol{w}}\left\{\frac{1}{n}\sum_{i=1}^n(\hat{y}_i (\boldsymbol{w})-y_i)^2\right\}}=0
\\

\Rightarrow{\frac{1}{n}\frac{\partial}{\partial\boldsymbol{w}}\left\{( \boldsymbol{Xw}-\boldsymbol{y})^2\right\}}=0
\\

\Rightarrow{\frac{1}{n}\frac{\partial}{\partial\boldsymbol{w}}\left\{ \boldsymbol{w}^{T}\boldsymbol{X^{T}Xw}-2\boldsymbol{w^{T}X^{T}y}+\boldsymbol{y^{T}y})\right\}}=0
\\

\Rightarrow2\boldsymbol{X^TXw}-2\boldsymbol{X^Ty}=0
\\

\Rightarrow\boldsymbol{X^TXw}=\boldsymbol{X^Ty}
\\

$Ax=b$の形、つまり$A=X^TX$と$b=X^Ty$を用いて$A$の逆行列を解いて$\boldsymbol{w}$を求める。式にすると

\boldsymbol{\hat{\boldsymbol{w}}}=(\boldsymbol{X^TX})^{-1}\boldsymbol{X}^{T}\boldsymbol{y}
\\

と書ける。

また、新しい入力値$X^{*}$に対する予測${\hat{y}}$は

\begin{align}
{\hat{y}}&=X^{*}\boldsymbol{\hat{w}}\\
&=X^{*}(\boldsymbol{X^TX})^{-1}\boldsymbol{X}^{T}\boldsymbol{y}
\end{align}

となる。

実装演習結果

m=1, n=10の時の例

import matplotlib.pyplot as plt
import numpy as np
from scipy import linalg

# an example with single parameter
# x: input  y:output

n=10;

x=np.arange(n);
y=np.random.rand(n);
y=x*y;

plt.plot(x,y,'.')

x=np.matrix(x);
x=x.T
x2=np.c_[np.ones(x.shape[0]),x]  # add ones to solve for w0

X=np.dot(x2.T,x2)      

y=np.matrix(y);
y=y.T
b=np.dot(x2.T,y)

# solve Xw=b

w=linalg.solve(X,b)    # w : weight
y2=np.dot(x2,w)        # y2: prediction

plt.plot(x,y)
plt.plot(x,y2)
plt.show()

参考文献

PythonでのMatrix計算においては、関連するウェブサイト等を参考にした。たとえば以下。
https://qiita.com/NNNiNiNNN/items/4fd5367f9ead6e5905a9

2.非線形回帰モデル

入力（説明変数）と出力（目的変数）との関係を多項式やガウス型基底関数などを用いた非線形回帰モデルで表現する。入力値$x_i$に対して非線形な基底関数$\boldsymbol{\phi(x_i)}$がk個あるモデルで計画行列はn x k次元となる。入力値に対して非線形であるが、係数$\boldsymbol{w}$については線形のモデルとすることで、係数$\boldsymbol{w}$の導出に際して線形回帰の場合と同様に扱うことができる。

\boldsymbol{\hat{w}}=(\boldsymbol{\Phi^{T}\Phi)^{-1}\Phi^{T}y}\\

{\hat{y}}=\Phi^{*}(\boldsymbol{\Phi^{T}\Phi)^{-1}\Phi^{T}y}

ただし、非線形モデルをより複雑化させてkを大きくしていってもデータとのフィッティングが向上していくとは限らない。一般にはシンプルな関係で表現できるものはシンプルなモデルを用いるのが良いとされる。
また個人的には、入力値のレンジが大きいデータに対して多項式モデルを採用すると、入力による基底関数の桁幅が大きくなり、解が不安定化した経験がある。また教師データの範囲外の入力値に対する不安定さ（外挿）は、多項式の項数が多いほど増していく。課題に応じた入力値のコンディショニングや多項式モデルの項数の検討などが必要と考えられる。

非線形回帰問題に限らないが、最適化関数に正則化項を加えることで、スムースな解や、スパースな解を求めるなどの工夫も一般的に用いられている。

教師データを用いた訓練誤差ではなく、教師データに用いないデータによる汎化誤差が最小化されているかを確認しながら基底関数の数、正則化項の選択やその項の重みを調整していく。

実装演習結果

３次式モデルでのフィッティング例を以下に示す。なお、教師データは非線形パターンを含むランダムなデータとした。

import matplotlib.pyplot as plt
import numpy as np
from scipy import linalg
n=50;

x=np.arange(n);
y=np.random.rand(n);
y=x*y;
yn=x*y+np.power(x,2)*y+np.power(x,3)*y

x=np.matrix(x);
x=x.T
x2=np.c_[np.ones(x.shape[0]),x,np.power(x,2),np.power(x,3)]  # nonlinear model

X=np.dot(x2.T,x2)      

yn=np.matrix(yn);
yn=yn.T
b=np.dot(x2.T,yn)

# solve Xw=b

w=linalg.solve(X,b)     # w : weight
yn2=np.dot(x2,w)        # y2: prediction

plt.plot(x,yn)
plt.plot(x,yn2)
plt.show()

3.ロジスティック回帰モデル

線形回帰と異なり目的変数が2値（0と1）などにする場合の教師あり問題に用いられる。
与えられた特徴量（説明変数）$\boldsymbol{x}$に対して目的変数（ラベル）yは0か1のどちらかとなるが、例えば1になる確率を$P(Y=1|X=\boldsymbol{x})$と表す。これは、分類問題における識別的アプローチ$p(C_k|x)$を直接モデル化する手法の1つであり、乗客の特徴に対する生死（1なら生存）や、表裏、成功不成功、犬または猫、などの問題がこれに当てはまる。

Y=0とY=1になる確率をまとめて表現するのにベルヌーイ分布が用いられる。

P(y)=p^y(1-p)^{1-y}

例えば、確率p=30％であたりの場合、y=1となるのは30％、y=0となるのは70％となる。
一方、当たり外れなどの2値データが複数あるときの尤度関数を求める際には、ベルヌーイ試行をn回行ったときの尤度関数（二項分布）で表される。$y_1, y_2, \cdots,y_n$が与えられたときに尤度関数が最大となるような確率pあるいはパラメータωを推定したい（最尤推定）。

P(y_1, y_2, \cdots,y_n;p)=\prod_{i=1}^{n}{p^{y_i}(1-p)^{1-y_i}}\\

ここで、確率pは入力値xと重みwの線形結合をシグモイド関数の入力としたもの$p=\sigma(\boldsymbol{\omega^Tx_i})$として表現できることから、上式は以下のようになる。

P(y_1, y_2, \cdots,y_n;\boldsymbol{\omega})=\prod_{i=1}^{n}{\sigma(\boldsymbol{\omega^Tx_i})^{y_i}(1-\sigma(\boldsymbol{\omega^Tx_i}))^{1-
y_i}}\\\\
=L(\boldsymbol{\omega})

最適化関数は、$L(\omega）$を最大にするパラメータ$\omega$を探す、つまりそれぞれの入力値の重みとシフト量を求めることになるが、$L(\omega）$はnが大きいと桁落ちしやすいことから$E(\omega)=-\log(L(\omega))$を最小化する問題に変形する。線形回帰の時と同様に、$E(\omega)$を$\omega$で微分して0になるときのパラメータ値を求める。逆行列を解析的に求めることは難しいため、大規模行列の逆問題を解く際などと同様に、ラインサーチ法（勾配法など）によって$E(\omega)$が小さくなる方向への繰り返し計算することで最適パラメータの組み合わせを求めていく。

\boldsymbol{\omega}^{(k+1)}=\boldsymbol{\omega}^{(k)}+\eta\sum_{i=1}^{n}{(y_i-p_i)\boldsymbol{x}_i}

ただし、$\eta$は学習率と呼ばれる。nが大きい大きなデータを扱う場合などは、すべてのデータに対して勾配計算をせずに、確率的に選択したデータに対して行うことで計算時間を小さくしたり容量不足を解決する。この手法は確率的勾配降下法と呼ばれている。

実装演習結果

ここでは、ロジスティック回帰の問題をつかむため、目的変数y=[1,0,0,1]、説明変数x=[[0,0,1,1],[0,1,0,1]]と[[1,0,0,1],[1,0,0,1]]の場合の分類精度を試した。
前者の場合は逆行列が求まらない（不能）となり分類できない結果となる。線形回帰問題でも同様であるが、x1とx2の数値が同じか違うかが目的変数を表していそうというデータや判断があれば、if$(x_{1i}=x_{2i})$などで$x_1$と$x_2$の一致不一致を表す新しい特徴量を入力とするなどの工夫が必要と考えられた。

import numpy as np
from sklearn.linear_model import LogisticRegression

x=np.array([[0,0,1,1],[0,1,0,1]])
y=np.array([1,0,0,1])
clf=LogisticRegression(max_iter=100,random_state=10)
clf.fit(x.T,y.T)
print(clf.score(x.T,y.T)) # Precision

x2=np.array([[1,0,0,1],[1,0,0,1]])
clf=LogisticRegression(max_iter=100,random_state=10)
clf.fit(x2.T,y.T)
print(clf.score(x2.T,y.T)) # Precision

参考文献

ロジスティック回帰の最適化関数の微分の定式化は例えば以下を参照
https://gihyo.jp/dev/serial/01/machine-learning/0019
勾配法などは、最適化問題の教科書を参照。
https://www.kyoritsu-pub.co.jp/kenpon/bookDetail/9784320017863

4.主成分分析

主成分分析は機械学習においては教師無し学習の1種であり次元圧縮の手法と説明されている。複数の特徴量を持つ元データ$\boldsymbol{x_i}$の特徴を最もよくとらえた低次元の射影を見つける。これは、次元圧縮後の分散を最大化することに相当する。

元データからその平均を引いたデータ行列を$X$とすると、その分散共分散行列は

Var(X)=\frac{1}{n}X^{T}X

線形変換後のベクトルを

\boldsymbol{s}_j=X\boldsymbol{a}_j

と置くと、線形変換後の分散は

Var(\boldsymbol{s}_j)=\frac{1}{n}\boldsymbol{s}_j^{T}\boldsymbol{s}_j\\
\\
=\boldsymbol{a}_j^{T}Var(X)\boldsymbol{a}_j

と表される。ノルム$a_j^{T}a_j=1$となる制約条件下で線形変換後の分散が最大となるのはラグランジュ未定乗数法を用い

E(\boldsymbol{a}_j)=\boldsymbol{a}_j^TVar(X)\boldsymbol{a}_j-\lambda(\boldsymbol{a}_j^T\boldsymbol{a}_j-1),\\

\frac{\partial E(\boldsymbol{a}_j)}{\partial \boldsymbol{a}_j}=0

から

Var(X)\boldsymbol{a}_j=\lambda\boldsymbol{a}_j

となる。射影先の第k主成分の分散$Var(\boldsymbol{s}_k)$は主成分に対応する固有値に一致する。この固有値は必ず0以上となり、固有ベクトルは互いに直行する。

Xを特異値分解し$X=U\Sigma V^T$とすると、

Var(X)＝\frac{1}{n}(V\Sigma V^{T})\\
        \\
　　　=V\left(
       \begin{array}{ccc}
          \sigma_1^2/n & & & \\
          & \sigma_2^2/n & & \\
          &  &  \ddots & \\
          &  &   &   \sigma_M^2/n 

       \end{array}
       \right) V^{T}

となる。これは、分散共分散行列の固有値が$\sigma_1^2/n, \sigma_2^2/n, \dots$であり、特異値分解することで射影先の分散、固有値が求まることを示す。

固有値(寄与率)を大きい順に並べ、第1～k主成分の分散の全分散に対する割合を累積寄与率とよび、k次元へ圧縮後に保存された情報の割合を表す。圧縮されたデータは分類や他の機械学習手法の入力に用いられる。

実装演習結果

ここでは以下のとおりscikit-learn並びにnumpyを用いて主成分分析を実施した。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from sklearn.datasets import load_iris

#load iris data

iris=load_iris()
data_iris=pd.DataFrame(iris.data, columns=iris.feature_names)
data_iris['target']=iris.target
data_iris=data_iris[0:100]

scaler=StandardScaler()
data_std=scaler.fit_transform(data_iris[iris.feature_names])
pca=PCA(n_components=4)
pca_transformed=pca.fit_transform(data_std)
plt.scatter(pca_transformed[:,0], pca_transformed[:,1],c=data_iris['target'])
print(pca_transformed.shape)
plt.plot(pca.explained_variance_ratio_)

dfc=data_iris

nDim=2                        # dimension after PCA
dfc=dfc-dfc.mean(axis=0)       
aa=np.cov(dfc, rowvar=False)   # covariance matrix
eig_val, eig_vec =np.linalg.eig(aa)  # eigen value and eigen vector
eigi=np.argsort(eig_val)[::-1]       # sort
eigv=v[:,eigi]
C=eigv[:,:nDim]
Out=np.dot(dfc,C)                # projection to reduced dimension

plt.scatter(Out[:,0],Out[:,1],c=data_iris['target'])

5．アルゴリズム

k 近傍法　

分類問題のための機械学習法であり、最近傍のデータをk個とり、それが最も多く所属するクラスに識別する手法である。kを大きくしていくことでより多くのデータによって所属が決定されるため、決定境界は滑らかになる。
以下を参照

https://dev.classmethod.jp/articles/2017ad_20171218_knn/
https://qiita.com/popotan322000/items/d4fbdfc02570673ac2c0

k 平均法　

与えられたデータをk個のクラスタに分類する、教師無し学習手法である。まず各クラスタ中心の初期値を設定し、各データから各中心との距離が最も近いクラスタを割り当てる。クラスタの平均ベクトルを計算して各中心の位置をアップデートし、繰り返し計算を行うものである。

実装演習結果

ここではscikit-learnのirisデータを用いたKNN(K-Nearest Neighbor Algorithm )のテストを実施した。

#KNN
import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
import matplotlib.pyplot as plt
#load iris data

iris=load_iris()
data_iris=pd.DataFrame(iris.data, columns=iris.feature_names)
data_iris['target']=iris.target
data_iris=data_iris[0:150]

Y = np.array(data_iris['target'])
X = np.array(data_iris[iris['feature_names']])

# divide data into learning and validating
from sklearn.model_selection import train_test_split
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2)

from sklearn.neighbors import KNeighborsClassifier

list_nn = []
list_score = []
for k in range(1, 51): 
  knc = KNeighborsClassifier(n_neighbors=k)
  knc.fit(X_train, Y_train)

  Y_pred = knc.predict(X_test)         # prediction
  score = knc.score(X_test, Y_test)    # evaluation
  print("[%d] score: {:.2f}".format(score) % k)

  list_nn.append(k)
  list_score.append(score)

plt.ylim(0.8, 1.0)
plt.xlabel("n_neighbors")
plt.ylabel("score")
plt.plot(list_nn, list_score)

6．サポートベクターマシン

一般に2クラス分類のために用いられる機械学習法である。線形判別関数と最も近いデータ点との距離（マージン）を最大とする線形判別関数を求める手法である。

n次元空間の線形式は

y=\omega_0+\boldsymbol{\omega^Tx}

とあわらせ、この符号が分類のラベルを示す。ラベル$y_i$の値は+1または-1とし正しく分類されているとすると、$y_i(\omega_0+\boldsymbol{\omega^Tx})>0$となる。そのとき、i番目の点から直線までの距離は以下となる。

\frac{y_i(\omega_0+\boldsymbol{\omega^Tx})}{||w||}

ここで、すべての点のなかで境界直線までの距離が一番近いものを最大化したいので、

Maximize \min_i\frac{y_i(\omega_0+\boldsymbol{\omega^Tx})}{||w||}

と表現される。これを変形していくと、$y_i(\omega_0+\boldsymbol{\omega^Tx})\geq1$の条件のもとで、$||w||^2$を最小化する問題に帰結する。ただし、この式だと例外値（正しく分類されない点）があると解くことができない。

こうした場合は、例外値もある程度対応するために、スラッグ変数$\xi$を導入し、例外値の許容度合を調整しながらマージンを最大化する。

以下参照
https://watlab-blog.com/2019/12/22/svm/
http://sudillap.hatenablog.com/entry/2013/04/08/235602

実装演習結果

ここではscikit-learnのirisデータを用いたSVMのテストを実施した。

## SVM
import numpy as np
import pandas as pd
from sklearn.datasets import load_iris
#load iris data

iris=load_iris()
data_iris=pd.DataFrame(iris.data, columns=iris.feature_names)
data_iris['target']=iris.target
data_iris=data_iris[0:150]

Y = np.array(data_iris['target'])
X = np.array(data_iris[iris['feature_names']])

# divide data into learning and validating
from sklearn.model_selection import train_test_split
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2)

from sklearn.svm import SVC
SVMmodel = SVC(gamma='scale')
SVMmodel.fit(X_train, Y_train)

Y_pred = SVMmodel.predict(X_test)

print(Y_test[:30])
print(Y_pred[:30])

参考文献

以下の記事が参考になった。
https://qiita.com/renesisu727/items/1cf367d17b46237f9c7a
また、加藤公一氏著、機械学習のエッセンスという本でも詳しく書かれている。
https://www.amazon.co.jp/exec/obidos/ASIN/4797393963/hatena-blog-22/