【ラビットチャレンジ レポート】機械学習編

#0. はじめに
本記事は、Study-AIが提供するラビット・チャレンジ（JDLA認定プログラム）のレポート記事である。今回は機械学習に関するレポートを記述する。

###参考文献
######線形回帰モデル
[1]Rで学ぶデータサイエンス 20（共立出版）

######主成分分析
[2]これなら分かる応用数学教室 最小二乗法からウェーブレットまで（金谷健一、共立出版）

######サポートベクターマシン
[3]はじめてのパターン認識（森北出版）

#1. 線形回帰モデル
###線形回帰とは？
ある入力（離散あるいは連続値）から出力（連続値）を予測する問題を回帰問題という。その中で、直線（超平面）で予測するモデルを、線形回帰モデルという。
一変数でモデルを予測することを単回帰、他変数でモデルを予測することを重回帰と呼ぶ。次には、重回帰を説明していく（単回帰は重回帰の$n=1$の場合なので、省略する）。

###重回帰
回帰分析は機械学習の基本とも言える手法なので、少し丁寧に記述する。
さて、${(x_{i1},x_{i2},\cdots,x_{iq},y_i)}\ 0\leq i\leq n$というデータに対して以下の値

\begin{align}
\mathrm{RSS}=\sum_{i=1}^n(y_i-a_0-\sum_{j=1}^qa_1x_{ij})^2=\sum_{i=1}^n\varepsilon_i^2
\end{align}

を最小にすることで${a_j}$の値を求めることを重回帰という。これにより以下の回帰式が得られる。

\begin{align}
y=\hat{a}_0+\sum_{j=1}^q\hat{a}_jx_{j}
\end{align}

これを用いると、各データは以下のように書ける。

\begin{align}
y_i=\hat{a}_0+\sum_{j=1}^q\hat{a}_jx_{j}+\varepsilon_i
\end{align}

ここで、${a_j}$のような値を回帰係数という。また、${\hat{a}_j}$におけるハットは、これらの値が推定値であることを表している。以下、この値を行列を用いて計算していく。まず、$n$個のデータを以下のように$\boldsymbol{X}$と$\boldsymbol{y}$で表す。

\begin{align}
\boldsymbol{X}=\begin{pmatrix}
1&x_{11}&x_{12}&\cdots&x_{1q}\\
1&x_{21}&x_{22}&\cdots&x_{2q}\\
\vdots&\vdots&\vdots&\ddots&\vdots\\
1&x_{n1}&x_{n2}&\cdots&x_{nq}\\
\end{pmatrix},\quad\boldsymbol{y}=\begin{pmatrix}
y_1\\
y_2\\
y_3\\
\vdots\\
y_n
\end{pmatrix}\end{align}

この$\boldsymbol{X}$は計画行列という。次に以下のようにおく。

\begin{align}
\boldsymbol{y}=\tilde{\boldsymbol{y}}+\boldsymbol{\varepsilon}
\end{align}

ただし

\begin{align}
\tilde{\boldsymbol{y}}=\begin{pmatrix}
a_0+\sum_{j=1}^qa_jx_{1j}\\
a_0+\sum_{j=1}^qa_jx_{2j}\\
a_0+\sum_{j=1}^qa_jx_{3j}\\
\vdots\\
a_0+\sum_{j=1}^qa_jx_{nj}\\
\end{pmatrix},\quad\boldsymbol{\varepsilon}=\begin{pmatrix}
\varepsilon_1\\
\varepsilon_2\\
\varepsilon_3\\
\vdots\\
\varepsilon_n
\end{pmatrix}
\end{align}

また、（回帰）係数ベクトルを

\begin{align}
\boldsymbol{a}=\begin{pmatrix}
a_0\\
a_1\\
\vdots\\
a_q
\end{pmatrix},\quad \hat{\boldsymbol{a}}=\begin{pmatrix}
\hat{a}_0\\
\hat{a}_1\\
\vdots\\
\hat{a}_q
\end{pmatrix}
\end{align}

とする。ここで、データ${(x_{i1},x_{i2},\cdots,x_{iq},y_i)}\ 0\leq i\leq n$に対する線形回帰モデルは行列表示で

\begin{align}
\boldsymbol{y}=\boldsymbol{Xa}+\boldsymbol{\varepsilon}
\end{align}

と書ける。ここで、誤差の2乗和:

\begin{align}
\Delta^2\equiv||\boldsymbol{\varepsilon}||^2=||\boldsymbol{y}-\boldsymbol{Xa}||^2=\sum_{i=1}^n(y_i-a_0-a_1x_{i1}-a_2x_{i2}-\cdots-a_qx_{iq})^2
\end{align}

を最小にするような$\boldsymbol{a}$の値を$\boldsymbol{a}$の推定量$\hat{\boldsymbol{a}}$として求めたい（最小二乗法）。最小値を与える解を解くために$\Delta^2$を$a_j(j=1,\cdots,q)$で偏微分すると

\begin{align}
\dfrac{\partial}{\partial \boldsymbol{a}}\Delta^2=\begin{pmatrix}
\dfrac{\partial\Delta^2}{\partial a_0}\\
\dfrac{\partial\Delta^2}{\partial a_1}\\
\vdots\\
\dfrac{\partial\Delta^2}{\partial a_q}\\
\end{pmatrix}=&\begin{pmatrix}
2\sum_{i=1}^n(y_i-a_0-a_1x_{i1}-\cdots-a_qx_{iq})(-1)\\
2\sum_{i=1}^n(y_i-a_0-a_1x_{i1}-\cdots-a_qx_{iq})(-x_{i1})\\
\vdots\\
2\sum_{i=1}^n(y_i-a_0-a_1x_{i1}-\cdots-a_qx_{iq})(-x_{iq})\\
\end{pmatrix}\\
=&-2\begin{pmatrix}
\sum_{i=1}^n1\cdot(y_i-a_0-a_1x_{i1}-\cdots-a_qx_{iq})\\
\sum_{i=1}^nx_{i1}\cdot(y_i-a_0-a_1x_{i1}-\cdots-a_qx_{iq})\\
\vdots\\
\sum_{i=1}^nx_{iq}\cdot(y_i-a_0-a_1x_{i1}-\cdots-a_qx_{iq})\\
\end{pmatrix}\\
\notag\\
=&-2\begin{pmatrix}
\sum_{i=1}^n1\cdot y_i-\sum_{i=1}^n1\cdot(a_0+a_1x_{i1}+\cdots+a_qx_{iq})\\
\sum_{i=1}^nx_{i1}y_i-\sum_{i=1}^nx_{i1}(a_0+a_1x_{i1}+\cdots+a_qx_{iq})\\
\vdots\\
\sum_{i=1}^nx_{iq}y_i-\sum_{i=1}^nx_{iq}(a_0+a_1x_{i1}+\cdots+a_qx_{iq})\\
\end{pmatrix}\\
\notag\\%
=&-2\begin{pmatrix}
\sum_{i=1}^n(\boldsymbol{X})_{i1}(\boldsymbol{y})_i-\sum_{i=1}^n(\boldsymbol{X})_{i1}(\boldsymbol{Xa})_i\\
\sum_{i=1}^n(\boldsymbol{X})_{i2}(\boldsymbol{y})_i-\sum_{i=1}^n(\boldsymbol{X})_{i2}(\boldsymbol{Xa})_i\\
\vdots\\
\sum_{i=1}^n(\boldsymbol{X})_{i(q+1)}(\boldsymbol{y})_i-\sum_{i=1}^n(\boldsymbol{X})_{i(q+1)}(\boldsymbol{Xa})_i\\
\end{pmatrix}\\
\notag\\%
=&-2\left\{\begin{pmatrix}
\sum_{i=1}^n(\boldsymbol{X})_{i1}(\boldsymbol{y})_i\\
\sum_{i=1}^n(\boldsymbol{X})_{i2}(\boldsymbol{y})_i\\
\vdots\\
\sum_{i=1}^n(\boldsymbol{X})_{i(q+1)}(\boldsymbol{y})_i
\end{pmatrix}
-\begin{pmatrix}
\sum_{i=1}^n(\boldsymbol{X})_{i1}(\boldsymbol{Xa})_i\\
\sum_{i=1}^n(\boldsymbol{X})_{i2}(\boldsymbol{Xa})_i\\
\vdots\\
\sum_{i=1}^n(\boldsymbol{X})_{i(q+1)}(\boldsymbol{Xa})_i
\end{pmatrix}\right\}\\
\notag\\%
=&-2\left\{\begin{pmatrix}
(\boldsymbol{X}^t\boldsymbol{y})_1\\
(\boldsymbol{X}^t\boldsymbol{y})_2\\
\vdots\\
(\boldsymbol{X}^t\boldsymbol{y})_{q+1}
\end{pmatrix}
-\begin{pmatrix}
(\boldsymbol{X}^t\boldsymbol{Xa})_1\\
(\boldsymbol{X}^t\boldsymbol{Xa})_2\\
\vdots\\
(\boldsymbol{X}^t\boldsymbol{Xa})_{q+1}
\end{pmatrix}\right\}\\
\notag\\%
=&-2(\boldsymbol{X}^t\boldsymbol{y}-\boldsymbol{X}^t\boldsymbol{Xa})
\end{align}

となる。これを最小にするには

\begin{align}
-2(\boldsymbol{X}^t\boldsymbol{y}-\boldsymbol{X}^t\boldsymbol{Xa})=\boldsymbol{O}\quad\Leftrightarrow\quad\boldsymbol{X}^t\boldsymbol{y}=\boldsymbol{X}^t\boldsymbol{Xa}
\end{align}

という式を解く必要がある（正規方程式）。ここで、$\boldsymbol{X}^t\boldsymbol{X}$が正則行列の場合を考えることによって、これを最小にする$\boldsymbol{a}$の推定値は

\begin{align}
\hat{\boldsymbol{a}}=(\boldsymbol{X}^t\boldsymbol{X})^{-1}\boldsymbol{X}^t\boldsymbol{y}
\end{align}

と書ける。

###実装演習結果
######データの準備
まずは今回使うデータフレームを読み込む

#ライブラリのインポート
from sklearn.datasets import load_boston
from pandas import DataFrame
import numpy as np

# ボストンデータを"boston"というインスタンスにインポート
boston = load_boston()

今回使うデータでは、'feature_names'にラベル名、'data'に説明変数、'target'に目的変数が格納されている。それぞれのデータの確認の仕方は以下のようにできる。

print(boston['feature_names']) #ラベルの確認
print(boston['data']) #各データの説明変数を確認
print(boston['target']) #各データの目的変数を確認

次にPandasのDataFｒameで、扱うデータをデータを読み込む。

#　説明変数らをDataFrameへ変換
df = DataFrame(data=boston.data, columns = boston.feature_names)

# 目的変数をDataFrameへ追加
df['PRICE'] = np.array(boston.target)

#データの中身を確認（最初の10行）
df.head(10)

最初の10行しか表示していないが、実際は506行データが存在する。

######単回帰分析
まずは必要なライブラリを読み込み、オブジェクトを生成する

# sklearnモジュールからLinearRegressionをインポート
from sklearn.linear_model import LinearRegression

# オブジェクト生成
model = LinearRegression()

例えばBostonデータの中で住宅価格'PRICE'（目的変数）に直結しそうな、住宅の部屋数'RM'で単回帰を行う。

# 説明変数dataを用意
data = df.loc[:, ['RM']].values

# 目的変数targetを用意
target = df.loc[:, 'PRICE'].values

そしてfit関数を用いて、dataでtargetを予測する。

# fit関数でパラメータ推定
model.fit(data, target)

では予測したパラメータを用いて、predict関数で部屋が$N$個の時にいくらなのか、予測していく。

#予測 -> 部屋数N個の場合はいくらか？？
print("部屋が1個の時" , model.predict([[1]]), "*1,000$")
print("部屋が2個の時" , model.predict([[2]]), "*1,000$")
print("部屋が3個の時" , model.predict([[3]]), "*1,000$")
print("部屋が4個の時" , model.predict([[4]]), "*1,000$")
print("部屋が5個の時" , model.predict([[5]]), "*1,000$")
print("部屋が6個の時" , model.predict([[6]]), "*1,000$")
print("部屋が7個の時" , model.predict([[7]]), "*1,000$")
print("部屋が8個の時" , model.predict([[8]]), "*1,000$")
print("部屋が9個の時" , model.predict([[9]]), "*1,000$")
print("部屋が10個の時" , model.predict([[10]]), "*1,000$")

部屋が1〜3この時、まともな予測ができていないことがわかる。この原因は、モデルの設計をミスしている（部屋数で家賃を表すのは適切ではなかった）、直線fittingがよくない、説明変数が足りていない、またはデータセットの中に1~3部屋のものがなく部屋数が少ない時のデータがうまく表現できなかった、など様々な要因が考えられる。

######重回帰分析
次に、重回帰分析を行っていく。データの準備はすでにできているので、新たなオブジェクトを生成するところから行う。今回は説明変数2つ（ボストン主要都市までの重み付き距離'DIS'と、部屋の数'RM'）で重回帰分析を行う。

# オブジェクトの生成
model3 = LinearRegression()

# データの準備
ｄata3 = df.loc[:, ['DIS', 'RM']].values
target3 = df.loc[:, 'PRICE'].values

# fit関数でパラメータ推定
model.fit(data3, target3)

予測は以下の通り。まずは'DIS'を固定で部屋数'RM'を変化させた結果から載せる。

print('距離が5、部屋数が1の時', model3.predict([[5.0, 1]]), '*1,000$')
print('距離が5、部屋数が2の時', model3.predict([[5.0, 2]]), '*1,000$')
print('距離が5、部屋数が3の時', model3.predict([[5.0, 3]]), '*1,000$')
print('距離が5、部屋数が4の時', model3.predict([[5.0, 4]]), '*1,000$')
print('距離が5、部屋数が5の時', model3.predict([[5.0, 5]]), '*1,000$')
print('距離が5、部屋数が6の時', model3.predict([[5.0, 6]]), '*1,000$')
print('距離が5、部屋数が7の時', model3.predict([[5.0, 7]]), '*1,000$')
print('距離が5、部屋数が8の時', model3.predict([[5.0, 8]]), '*1,000$')
print('距離が5、部屋数が9の時', model3.predict([[5.0, 9]]), '*1,000$')
print('距離が5、部屋数が10の時', model3.predict([[5.0, 10]]), '*1,000$')

予想通り、部屋数が増えると価格も増加した。しかし、相変わらず部屋数が1〜3の時はマイナスのままであった。
次に'RM'を固定で、'DIS'を変化させた結果を載せる。

print('距離が1、部屋数が5の時', model3.predict([[1.0, 5]]), '*1,000$')
print('距離が2、部屋数が5の時', model3.predict([[2.0, 5]]), '*1,000$')
print('距離が3、部屋数が5の時', model3.predict([[3.0, 5]]), '*1,000$')
print('距離が4、部屋数が5の時', model3.predict([[4.0, 5]]), '*1,000$')
print('距離が5、部屋数が5の時', model3.predict([[5.0, 5]]), '*1,000$')
print('距離が6、部屋数が5の時', model3.predict([[6.0, 5]]), '*1,000$')
print('距離が7、部屋数が5の時', model3.predict([[7.0, 5]]), '*1,000$')
print('距離が8、部屋数が5の時', model3.predict([[8.0, 5]]), '*1,000$')
print('距離が9、部屋数が5の時', model3.predict([[9.0, 5]]), '*1,000$')
print('距離が10、部屋数が5の時', model3.predict([[10.0, 5]]), '*1,000$')

こちらは予想に反して、重み付き距離が増加すると、家賃も増加した。

#2. 非線形回帰モデル
1章で説明した線形回帰モデルは、回帰式が

\begin{align}
y=\hat{a}_0+\sum_{j=1}^q\hat{a}_jx_{j}
\end{align}

と、$\sum_{j=1}^q\hat{a}_jx_j$の部分が線形で表されているので、線形で表現できるデータしか予測できないという難点があった。そこで、線形の部分を任意の関数に変えることで、

y=\hat{a}_0+\sum_{j=1}^q\hat{a}_j\phi_ｊ(x)

と、非線形の傾向を持つデータも予測できるようになる。導出のエッセンスは、線形でも非線形でも変わらず、誤差関数を最小にするような推定値$\hat{a}_j$を求めれば良い。例えば非線形関数に使われる関数として、多項式関数

\phi_j(x)=x^j

や、ガウス型基底関数

\phi_j(x)=\exp\left\{-\dfrac{(x-\mu_j)^2}{\sigma^2}\right\}

がある。

###実装演習結果
######データの準備
まずはライブラリを読み込み、seabornで作図の準備をする。

import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

#seaborn設定
sns.set()
#背景変更
sns.set_style("darkgrid", {'grid.linestyle': '--'})
#大きさ(スケール変更)
sns.set_context("paper")

次に、データを準備する。

n=100

def true_func(x):
    z = 1-48*x+218*x**2-315*x**3+145*x**4
    return z 

def linear_func(x):
    z = x
    return z 

# 真の関数からノイズを伴うデータを生成

# 真の関数からデータ生成
data = np.random.rand(n).astype(np.float32)
data = np.sort(data)
target = true_func(data)

# 　ノイズを加える
noise = 0.5 * np.random.randn(n) 
target = target  + noise

# ノイズ付きデータを描画

plt.scatter(data, target)

plt.title('NonLinear Regression')
plt.legend(loc=2)

以下のようなデータを用いる。4次関数（$145x^4-315x^3+218x^2-48x+1$）を計算する関数True_func(x)からノイズを発生させたものを用いる。

######非線形なデータを線形回帰で予測する。
まず、1章で行った線形回帰を行ってみる。

from sklearn.linear_model import LinearRegression

clf = LinearRegression()
data = data.reshape(-1,1)
target = target.reshape(-1,1)
clf.fit(data, target)

p_lin = clf.predict(data)

plt.scatter(data, target, label='data')
plt.plot(data, p_lin, color='darkorange', marker='', linestyle='-', linewidth=1, markersize=6, label='linear regression')
plt.legend()
print(clf.score(data, target))

結果は以下の通りである。

まったく予測ができていない。トレンドとして徐々に右肩上がりになっていることはわかるが、予測としては物足りない。

######非線形データを非線形回帰で予測する
そこで、非線形関数をでフィッティングしてみる。

from sklearn.kernel_ridge import KernelRidge

clf = KernelRidge(alpha=0.0002, kernel='rbf')
clf.fit(data, target)

p_kridge = clf.predict(data)

plt.scatter(data, target, color='blue', label='data')

plt.plot(data, p_kridge, color='orange', linestyle='-', linewidth=3, markersize=6, label='kernel ridge')
plt.legend()
#plt.plot(data, p, color='orange', marker='o', linestyle='-', linewidth=1, markersize=6)

こちらはノイズが乗って少しバラバラになったデータの傾向をしっかり捉え、表現できている。このように非線形回帰でしか予測できないデータもあるため、扱うデータによってどの回帰方法が適しているのかをしっかりと考えて用いなければならない。

#3. ロジスティック回帰モデル
ロジスティック回帰は、分類問題を解くのに適したモデルである。ロジスティック回帰は、ある入力（説明変数または特徴量の集まり）

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

からクラスに分類するための出力（目的変数）

y\in\{0, 1\}

を得るモデルである。求めるパラメータを

\boldsymbol{\omega}=(\omega_1, \omega_2, \cdots, \omega_m)^T\in\mathbb{R}^m, \text{and }\omega_0

とすると、出力は以下のように書ける。

\sigma(\hat{y})=\sigma(\boldsymbol{\omega}^T \boldsymbol{x}+\omega_0)

ここで、$\sigma(x)$はシグモイド関数である。

######シグモイド関数
入力に実数をとり、出力は必ず0〜1の値をとる、単調増加関数である。

\sigma(x)=\dfrac{1}{1+\exp(-ax)}

パラメータ$a$をずっと増加させると、ステップ関数に近づく。シグモイド関数の大切な性質として、微分しても自分自身で表現することが可能である。また、微分すると、

\dfrac{\partial \sigma(x)}{\partial x}=\sigma(x)（1-\sigma(x)）\leq\dfrac{1}{4}

となる。この性質により、深層学習する際に勾配消失問題を引き起こす。

######推定値の決定
ロジスティック回帰のパラメータの推定値は、尤度関数

\begin{align}
L(\boldsymbol{\omega})=&P(y_1, y_2, \cdots, y_n|\omega_0, \omega_1, \cdots, \omega_m)\\=&\prod^n_{i=1}p_i^{y_i}(1-p_i)^{1-y_i}\\
=&\prod^n_{i=1}\sigma(\boldsymbol{\omega}^T\boldsymbol{x_i})^{y_i}(1-\boldsymbol{\omega}^T\boldsymbol{x_i})^{1-y_i}
\end{align}

を最大化するパラメータを探索する。しかし実際は、尤度関数の負の対数尤度（negative log-likelyhood）$-\log L(\boldsymbol{\omega})$を最小化する。対数を取る理由の一つに実装上の問題があり、$0\leq p\leq1$を何度も掛けると、間違いなく桁落ちしてしまうのを回避するためである。

###実装演習結果
######データの準備
まずはライブラリの準備を行う。

import pandas as pd
from pandas import DataFrame
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from sklearn.linear_model import LogisticRegression

今回使うデータはタイタニックデータと呼ばれ、タイタニック号に乗った人に関する特徴量（年齢や客室のグレード等）と、生存できたかできなかったかという目的関数が記録されている。このデータから、ある特徴量を持つ人は生存できるかできないかを予測する（生存する確率を求める）ことができるモデルを作成するのが我々の目的である。
今回は事前にダウンロードしたタイタニックデータを読み込んで使うため、以下のコマンドを用いて格納場所から読み込む。

titanic_df = pd.read_csv('[格納場所のパス]/titanic_train.csv')
titanic_df.head(10)

このようなデータになっている。

この中の「Survived」が1ならこの人は生存した、0なら残念ながら生存することができなかったことを表す。

######データの前処理
さて、今回用いるデータは以下の点でそのままでは使えない。

1. データに文字列が含まれる
2. データにnullが含まれる
   まずは1.から解決していく。今回は必要のないデータを削除し、必要なデータ（例えば性別）を数値に変換することを行う。

#不要なデータの削除
titanic_df.drop(['PassengerId', 'Name', 'Ticket', 'Cabin'], axis=1, inplace=True)
#性別データを数値（女性なら0、男性なら1）に変換
titanic_df['Gender'] = titanic_df['Sex'].map({'female': 0, 'male': 1}).astype(int)

次に、2.を解決する。nullが含まれるデータの解決法は色々ある。例えばnullが含まれるデータを使わない、nullを適当な数値で埋める、nullをありそうな数値で埋める（補間する）というものである。今回はfillna()関数を用いて、年齢がnullの場合、データ全体の年齢の平均で埋めるという方法をとる。

#Ageカラムのnullを平均で補完（Ageに直接上書きせずに、新たな列AgeFillを作る）
titanic_df['AgeFill'] = titanic_df['Age'].fillna(titanic_df['Age'].mean())

次から、今回前処理したデータを用いて予測をしていく。

######１変数のロジスティック回帰
ここでは、チケットの価格「Fare」から生死を予測していく。高い部屋にいるほうが設備がちゃんとしていて生存しやすいだろうという直感的な予測とあっているのかどうかを検証する。

data1 = titanic_df.loc[:, ["Fare"]].values        # 運賃のリストの作成
label1 =  titanic_df.loc[:,["Survived"]].values   # 生死のリストを作成

そしてロジスティック回帰のオブジェクトを作成し、上記データを入れ予測を開始する。ここでは最尤推定により推定値を決定している。

model=LogisticRegression()
model.fit(data1, label1)

では上記で予測したモデルを用いて、ある運賃の時にその人は生存できるのかどうかを試してみる。

for i in [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300]:
  a = model.predict([[i]])
  if a == 0:
    print('運賃が{}の時、その部屋に泊まっている乗客は死亡する'.format(i))
  else:
    print('運賃が{}の時、その部屋に泊まっている乗客は生存する'.format(i))

このような結果になった。やはり高い部屋に泊まっている乗客は生存しやすいということがわかった。ロジスティック回帰の、他の分類手法と異なる面白いところは、そのクラスに分類される確率を求められる点である。例えば今回の場合、生存する確率を検証すると、

for i in [10, 20, 30, 40, 50, 60, 70, 80, 90, 100, 200, 300]:
  a = model.predict_proba([[i]])
  print('運賃が{}の時、その部屋に泊まっている乗客が生存する確率は{:.3f}'.format(i, a[0][1]))

上記のようになっていることがわかる。

######2変数のロジスティック回帰
今度は2変数として、「Pclass_Gender」と年齢で用いる。「Pclass_Gender」とは、Pclass（乗客の社会階級）とGender（年齢）を以下のように組み合わせる。

titanic_df['Pclass_Gender'] = titanic_df['Pclass'] + titanic_df['Gender']

これは、上流階級であるほど数字が小さくなるのに加え、男性は女性より値が1大きくなる。上流階級ほど生存しやすく、さらに女性の方が生存しやすいということを表した新たな変数である。さて、これらの変数を用いると、以下のような結果が得られる。

np.random.seed = 0

xmin, xmax = -5, 85
ymin, ymax = 0.5, 4.5

index_survived = titanic_df[titanic_df["Survived"]==0].index
index_notsurvived = titanic_df[titanic_df["Survived"]==1].index

from matplotlib.colors import ListedColormap
fig, ax = plt.subplots()
cm = plt.cm.RdBu
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
sc = ax.scatter(titanic_df.loc[index_survived, 'AgeFill'],
                titanic_df.loc[index_survived, 'Pclass_Gender']+(np.random.rand(len(index_survived))-0.5)*0.1,
                color='r', label='Not Survived', alpha=0.3)
sc = ax.scatter(titanic_df.loc[index_notsurvived, 'AgeFill'],
                titanic_df.loc[index_notsurvived, 'Pclass_Gender']+(np.random.rand(len(index_notsurvived))-0.5)*0.1,
                color='b', label='Survived', alpha=0.3)
ax.set_xlabel('AgeFill')
ax.set_ylabel('Pclass_Gender')
ax.set_xlim(xmin, xmax)
ax.set_ylim(ymin, ymax)
ax.legend(bbox_to_anchor=(1.4, 1.03))

この結果を見ると、Pclass_Genderが小さい（上流階級）ほど生存しやすく、大きい（下流階級）ほど死亡しやすくなる。また、男性よりも女性の方が生存しやすい事がわかる。年齢に関しては、歳をとるほど死亡しやすく、若いほど生存しやすい傾向にあるようだ。

#4. 主成分分析
多数の多次元データを解析して意味のある量を引き出す手法を主成分分析という。主に次元削除などに使われる。主成分分析は、データの分散が最大となる方向の軸を取る。ではどうすれば、多次元で分散が最大になる方向を知ることができるのか。

######分散・共分散行列を対角化する
特徴量$i$に対する分散を

\begin{align}
V_{ii}=Var\big(f(X=x)\big)=E\Big[\Big(f(x)-E\big(f(x)\big)\Big)^2\Big]=E\big(f(x)^2\big)-\Big(E\big(f(x)\big)\Big)^2
\end{align}

とし、2つの特徴量$i$と$j$間の共分散を、

\begin{align}
C_{ij}=Cov\big(f(X=x), g(Y=y)\big)&=E\Big[\Big(f(x)-E\big(f(x)\big)\Big)\Big(g(y)-E\big(g(y)\big)\Big)\Big]\\
&=E\big(f(x)g(x)\big)-E\big(f(x)\big)E\big(f(y)\big)
\end{align}

とする。そのとき、分散・共分散行列は次のように表される。

\begin{align}
\boldsymbol{V}=\begin{pmatrix}
V_{11}&C_{12}&\cdots&C_{1N}\\
C_{21}&V_{22}&\cdots&C_{2N}\\
\vdots&\cdots&\ddots&\vdots\\
C_{N1}&C_{N2}&\cdots&V_{NN}
\end{pmatrix}
\end{align}

この分散・共分散行列$\boldsymbol{V}$を対角化すると、固有値$\lambda_1, \lambda_2, \cdots, \lambda_N$とそれに対応する固有ベクトル$\boldsymbol{v}_1, \boldsymbol{v}_2, \cdots, \boldsymbol{v}_N$が得られる。実はこの中の最大固有値（ここでは仮に$\lambda_1$とする）に対応する固有ベクトル（ここでは$\boldsymbol{v}_1$）の軸が、分散が最大になる軸であり、主成分分析では第1主成分となる。もし固有値が降順に並んでいたとしたら、第$k$主成分に対応する軸は、$\boldsymbol{v}_k$で表される。

###実装演習結果
######データの準備

まずはライブラリを読み込む。

import pandas as pd
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegressionCV
from sklearn.metrics import confusion_matrix
from sklearn.decomposition import PCA
import matplotlib.pyplot as plt

次に、データを読み込む。今回はローカルにあるデータを読み込んでいく。

cancer_df = pd.read_csv('[データを格納しているフォルダのパス]/data/cancer.csv')

今回はデータが多いためデータの提示は割愛するが、33次元、569データあるものを用いる。

# 目的変数の抽出（分類のための変数diagnosisがMなら1、それ以外（B）なら0とする）
y = cancer_df.diagnosis.apply(lambda d: 1 if d == 'M' else 0)

# 説明変数の抽出
X = cancer_df.loc[:, 'radius_mean':]

# 学習用とテスト用でデータを分離
X_train, X_test, y_train, y_test = train_test_split(X, y, random_state=0)

# 標準化
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled = scaler.transform(X_test)

######主成分分析
主成分分析を行う。オブジェクトを生成し、結果を予測(fit）する。

#どれくらい主成分が情報量を持つか？
#第一主成分は43％、第二主成分は19%程度なので、2次元でも62%ほど説明できる。
pca = PCA(n_components=30)
pca.fit(X_train_scaled)
plt.bar([n for n in range(1, len(pca.explained_variance_ratio_)+1)], pca.explained_variance_ratio_)

pca.explained_variance_ratio_は固有値を規格化し、第$i$主成分がどれくらいの説明力を持つかを示す。結果は以下のようになる。

横軸は第$i$主成分、縦軸が比率を表す。これを見ると、第1主成分は43％、第2主成分は19%程度なので、2次元でも62%ほど説明できる。第6主成分まで使うと約90%を説明できるが、説明変数が多くなると軸の理解が難しくなったり、図示による解釈ができなくなるので今回は2次元に抑えておく。
では、実際に第2主成分までを用いて、分類の様子を見てみる。

# 次元数2まで圧縮
pca = PCA(n_components=2)
X_train_pca = pca.fit_transform(X_train_scaled)
print('X_train_pca shape: {}'.format(X_train_pca.shape))
# -> X_train_pca shape: (426, 2)

# 寄与率
print('explained variance ratio: {}'.format(pca.explained_variance_ratio_))
# -> explained variance ratio: [ 0.43315126  0.19586506]

# 散布図にプロット
temp = pd.DataFrame(X_train_pca)
temp['Outcome'] = y_train.values
b = temp[temp['Outcome'] == 0]
m = temp[temp['Outcome'] == 1]
plt.scatter(x=b[0], y=b[1], marker='o') # 良性は○でマーク
plt.scatter(x=m[0], y=m[1], marker='^') # 悪性は△でマーク
plt.xlabel('PC 1') # 第1主成分をx軸
plt.ylabel('PC 2') # 第2主成分をy軸

このように決定境界が曖昧になってしまうが、次元削除をすると説明がしやすくなることがわかる。通常はPC1軸の成分はどの説明変数が多く寄与していて、どんな特徴量を表すかなどの解析も行う。

#5. kクラスタリング手法
###k近傍法（kNN）
新しいデータ点を分類するときに、そのデータからもっとも近くにあるデータ（最近傍データ）を$k$個取ってきて、それらがもっとも多く所属するクラスに識別する手法である。この手法では、$k$がハイパーパラメータであり、$k$を大きくすると決定境界は滑らかになる（この結果はこの後のハンズオンで示す）。

###実装演習結果（kNN）
まずは必要なライブラリをimportする。

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

今回はサンプルデータは乱数を用いて作成するため、サンプルデータを作成する関数 gen_data() を実装する。

def gen_data():
    x0 = np.random.normal(size=50).reshape(-1, 2) - 1
    x1 = np.random.normal(size=50).reshape(-1, 2) + 1.
    x_train = np.concatenate([x0, x1])
    y_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
    return x_train, y_train

# サンプルデータの作成
X_train, ys_train = gen_data()
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)

今回は上記のデータになった（seedを固定していない乱数を使っているので、毎回生成されるデータが違う）。続いて、このサンプルデータの中に新たなデータを置いたらどっちのクラスに割り振られるかを予測していく。予測するデータ点との、距離が最も近い$k$個の、訓練データのラベルの最頻値を割り当てる関数を作成していく。

def distance(x1, x2):
    return np.sum((x1 - x2)**2, axis=1)

def knc_predict(n_neighbors, x_train, y_train, X_test):
    y_pred = np.empty(len(X_test), dtype=y_train.dtype)
    for i, x in enumerate(X_test):
        distances = distance(x, X_train)
        nearest_index = distances.argsort()[:n_neighbors]
        mode, _ = stats.mode(y_train[nearest_index])
        y_pred[i] = mode
    return y_pred

def plt_resut(x_train, y_train, y_pred):
    xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
    xx = np.array([xx0, xx1]).reshape(2, -1).T
    plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
    plt.contourf(xx0, xx1, y_pred.reshape(100, 100).astype(dtype=np.float), alpha=0.2, levels=np.linspace(0, 1, 3))

これらの関数を用いて、決定境界を作成する。

n_neighbors = 1 # <- ここがkNNのkである。

xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
X_test = np.array([xx0, xx1]).reshape(2, -1).T

y_pred = knc_predict(n_neighbors, X_train, ys_train, X_test)
plt_resut(X_train, ys_train, y_pred)

kNNの$k$の値を変化させると、以下のようになる。

データが密集している部分をよく見てみると、決定境界が滑らかになっていることがわかる。$k=1$の時は、1つデータが黄色いクラスに飛び出ているデータ（外れ値）に決定境界が引っ張られている。しかしそれが$k=3$、$k=10$と増やしていくと、決定境界が外れ値を無視するようになる。kNN法はパラメータ$k$が命であることがよくわかった。

###k-means法（k-平均法）
k-means法は教師なし学習のクラスタリング手法の1種である。与えられたデータを$k$個のクラスタに分類する。$k$は自分で決めるハイパーパラメータである。手順はシンプルであり、以下のアルゴリズムである。

1. 各クラスタの中心となる初期座標を設定する
2. 各データ点に対して、各クラスタの中心との距離を計算し、もっとも距離が近いデータ点をそのクラスタに所属するデータとして割り当てる
3. 各クラスタの平均ベクトル（中心）を計算する
4. 各クラスタの中心が収束するまで2. と3. の処理を繰り返す

このようにしてクラスタリングしていく。k-means法において重要なのは、クラスタをいくつにするかのハイパーパラメータ$k$の決定と、上記手順1. の初期座標の設定である。この2つで全てが決まる。

###実装演習結果（k-means法）
######データの準備
まずは必要なライブラリを読み込む。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn import cluster, preprocessing, datasets
from sklearn.cluster import KMeans

次に、今回用いるデータセット（wine）を読み込む。

wine = datasets.load_wine()
X = wine.data
y = wine.target

X.shape # -> (178, 13)

Xにはワインの特徴を表すデータを代入した。特徴量が13のデータが178個あることがわかる。yには目的変数（ワインのクラスを表す0〜2の数字）を入れる。

######k-meansの実行
KMeansオブジェクトを準備する。

#n_clustersはクラスター数を表す
model = KMeans(n_clusters=3)
labels = model.fit_predict(X)
df = pd.DataFrame({'labels': labels})

# 結果を見るために、クラスにいくつのデータが属しているかを調べる
def species_label(theta):
    if theta == 0:
        return wine.target_names[0]
    if theta == 1:
        return wine.target_names[1]
    if theta == 2:
        return wine.target_names[2]

df['species'] = [species_label(theta) for theta in wine.target]
pd.crosstab(df['labels'], df['species'])

「KMeans(n_clusters=3)」のn_clustersは、いくつのクラスターに分類したいか決めるハイパーパラメータである。今回は3クラスしかないからn_clusters=3としているが、未知のデータではこの値を慎重に決めなければならない。
結果は以下のようになった。

species	class_0	class_1	class_2
labels
0	13	20	29
1	46	1	0
2	0	50	19
全くうまく分類できていないことがわかる。特徴量が不十分であるのか、余計な特徴量が混ざってしまいクラスターを形成できてない可能性があることがわかる。

#6. サポートベクターマシン
###概要
サポートベクターマシン（SVM）とは、教師あり学習において回帰・分類。外れ値の検出を行う方法の一つである。SVMでは例えば2つの分類に対して、データ点をもっとも引き離すような境界線を引くことが目的である。特徴量（説明変数）が2つなら直線、3つなら平面が境界となるが、4つ以上になると超平面で表されるため図示はできない。

######サポートベクター（サポートベクトル）
サポートベクトルとは境界線にもっとも近いデータ点のことである。新しいデータが入力された時の誤判定を防ぐため、境界に近いデータであっても境界からできるだけ離すことが重要になる。サポートベクトルと境界との距離をマージンといい、SVMではこのマージンを最大化する。

###導出
######目的
最大マージン$D_\max$を実現する2クラス線形識別関数の学習法がSVMである。ここでは、マージン$D$が最大マージン（$D_\max$）となるような$\boldsymbol{w}$を導出する。

######導出
標準座標系を考える。クラスラベル付き学習データの集合を$\mathbb{D}=\{(ｔ_i, x_i)\} (i=1, 2, \cdots, N)$とする。$t_i=\{-1, +1\}$は教師データであり、学習データ$\boldsymbol{x}_i\in\mathbb{R}^d$がどちらのクラスに属するのかを指定する。
以下の図のように、線形識別関数のマージンを$\kappa$とすれば、全ての学習データで

|\boldsymbol{w}^T\boldsymbol{x}_i+b|\ge\kappa

が成り立つ。

係数ベクトルとバイアス項をマージン$\kappa$で正規化したものを改めて$\boldsymbol{w}$と$b$とすると、線形識別関数は、

\left\{
\begin{array}{ll}
\boldsymbol{w}^T\boldsymbol{x}_i+b\ge+1 & (t_i=+1) \\
\boldsymbol{w}^T\boldsymbol{x}_i+b\le-1 & (t_i=-1)
\end{array}
\right.

となる。この場合分けは、

t_i(\boldsymbol{w}^T\boldsymbol{x}_i+b)\ge1

のようにまとめることができる。
クラス間マージンは、各クラスのデータを$\boldsymbol{w}$の方向へ射影した長さの差の最小値

\rho(\boldsymbol{w}, b)=\min_{\boldsymbol{x}_i\in C_{t_i=+1}}\dfrac{\boldsymbol{w}^T\boldsymbol{x}}{||\boldsymbol{w}||}
-\max_{\boldsymbol{x}_i\in C_{t_i=-1}}\dfrac{\boldsymbol{w}^T\boldsymbol{x}}{||\boldsymbol{w}||}=\dfrac{1-b}{||\boldsymbol{w}||}-\dfrac{-1-b}{||\boldsymbol{w}||}=\dfrac{2}{||\boldsymbol{w}||}

で与えられる。最適な超平面の式を$\boldsymbol{w}_0^T\boldsymbol{x}_i+b_0=0$とすれば、この超平面は最大クラスマージン

\rho(\boldsymbol{w}_0, b_0)=\max_\boldsymbol{w}\rho(\boldsymbol{w}, b)

を与える。最大マージン$D_\max$は最大クラス間マージンの$1/2$で与えられる。したがって、最適識別超平面は$t_i(\boldsymbol{w}^T\boldsymbol{x}_i+b)\ge1 (i=1, 2 \cdots, N)$の制約下で、$\boldsymbol{w}$のノルムを最小にする解

\boldsymbol{w}_0=\min||\boldsymbol{w}||

として求めることができる。

###いろいろなSVMの手法
SVMでは、直線・平面・超平面を境界としてデータを分離するが、このことを線形分離という。しかし実際のデータにおいては、線形分離が可能なケースは多くない（何らかのノイズによって境界が不明瞭であったり、そもそもデータが線形分離可能な形をしていない場合があるため）。このような場合に使う手法に、ソフトマージンSVMとカーネル法がある。

######ソフトマージンSVM
ソフトマージンSVMは、誤差を認めるとより良い線形分離ができるようなデータに用いる。通常のSVM（ハードマージンSVMという）とは異なり、マージンにデータが入ることを許すのが特徴である。ただし、データがマージンに入った場合はペナルティを与え、マージンの最大化とペナルティの最小化を同時に行うことで、できるだけうまく分離するような境界を見つける。データには通常ノイズが含まれるため、実際の機械学習では主にソフトマージンSVMが用いられる。

######カーネル法
カーネル法は、クラスの境界が線形では表せない（曲線を使うしかない）場合に用いる。クラスの境界が曲面場合は、データから新しい特徴量を作って、うまく線形分離可能になるようにプロットする。この作業を、**「線形分離可能な高次元特徴空間に写像する」**という。データに存在する特徴量から、うまく線形分離可能になるような特徴量を新たに作り出す手法として、カーネル関数を用いる。カーネル関数としてはガウシアン（RBF）カーネルや、多項式カーネルなどがある。カーネル関数を使いつつ計算量を減らす工夫は、カーネルトリックと呼ばれる。

###実装演習結果
######通常のSVMのデータの準備
ライブラリをインポートし、今回用いるデータを正規分布から生成される乱数により作成する。

%matplotlib inline
import numpy as np
import matplotlib.pyplot as plt

def gen_data():
    x0 = np.random.normal(size=50).reshape(-1, 2) - 2.
    x1 = np.random.normal(size=50).reshape(-1, 2) + 2.
    X_train = np.concatenate([x0, x1])
    ys_train = np.concatenate([np.zeros(25), np.ones(25)]).astype(np.int)
    return X_train, ys_train

X_train, ys_train = gen_data()
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)

今回はこのようなデータになった。SVMにより、このデータに境界線を引くことを目的にする。

######通常のSVMの実行
下記のように境界線を引くことができる。目的は、マージンの最大化である。

t = np.where(ys_train == 1.0, 1.0, -1.0)

n_samples = len(X_train)
# 線形カーネル
K = X_train.dot(X_train.T)

eta1 = 0.01
eta2 = 0.001
n_iter = 500

H = np.outer(t, t) * K

a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.where(a > 0, a, 0)


index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()


xx0, xx1 = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * sv.dot(X_test[i])
y_pred = np.sign(y_project)


# 訓練データを可視化
plt.scatter(X_train[:, 0], X_train[:, 1], c=ys_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
# 領域を可視化
#plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])

下記のような結果になった。

本データは2クラスのデータがはっきり分かれていたのでハードマージンSVMでの対応がベストであったが、境界が曖昧だった場合はソフトマージンSVM、直線では境界線が引けないのであればカーネル法などを用いて対応する必要がある。

######線形分離不可なデータの準備
線形分離不可な、曲線でないと境界線が引けないようなデータを生成する。

factor = .2
n_samples = 50
linspace = np.linspace(0, 2 * np.pi, n_samples // 2 + 1)[:-1]
outer_circ_x = np.cos(linspace)
outer_circ_y = np.sin(linspace)
inner_circ_x = outer_circ_x * factor
inner_circ_y = outer_circ_y * factor

X = np.vstack((np.append(outer_circ_x, inner_circ_x),
               np.append(outer_circ_y, inner_circ_y))).T
y = np.hstack([np.zeros(n_samples // 2, dtype=np.intp),
               np.ones(n_samples // 2, dtype=np.intp)])
X += np.random.normal(scale=0.15, size=X.shape)
x_train = X
y_train = y

plt.scatter(x_train[:,0], x_train[:,1], c=y_train)

今回はこのようなデータになった。

######線形分離不可な場合のSVMの実行（カーネル法）
では、実際にRBFカーネル（ガウシアンカーネル）法を用いて高次元に飛ばし、境界線を引いてみる。

def rbf(u, v):
        sigma = 0.8
        return np.exp(-0.5 * ((u - v)**2).sum() / sigma**2)
    
X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)

n_samples = len(X_train)
# RBFカーネル
K = np.zeros((n_samples, n_samples))
for i in range(n_samples):
    for j in range(n_samples):
        K[i, j] = rbf(X_train[i], X_train[j])

eta1 = 0.01
eta2 = 0.001
n_iter = 5000

H = np.outer(t, t) * K

a = np.ones(n_samples)
for _ in range(n_iter):
    grad = 1 - H.dot(a)
    a += eta1 * grad
    a -= eta2 * a.dot(t) * t
    a = np.where(a > 0, a, 0)


index = a > 1e-6
support_vectors = X_train[index]
support_vector_t = t[index]
support_vector_a = a[index]

term2 = K[index][:, index].dot(support_vector_a * support_vector_t)
b = (support_vector_t - term2).mean()


xx0, xx1 = np.meshgrid(np.linspace(-1.5, 1.5, 100), np.linspace(-1.5, 1.5, 100))
xx = np.array([xx0, xx1]).reshape(2, -1).T

X_test = xx
y_project = np.ones(len(X_test)) * b
for i in range(len(X_test)):
    for a, sv_t, sv in zip(support_vector_a, support_vector_t, support_vectors):
        y_project[i] += a * sv_t * rbf(X_test[i], sv)
y_pred = np.sign(y_project)


# 訓練データを可視化
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
# サポートベクトルを可視化
plt.scatter(support_vectors[:, 0], support_vectors[:, 1],
                    s=100, facecolors='none', edgecolors='k')
# 領域を可視化
plt.contourf(xx0, xx1, y_pred.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
# マージンと決定境界を可視化
plt.contour(xx0, xx1, y_project.reshape(100, 100), colors='k',
                     levels=[-1, 0, 1], alpha=0.5, linestyles=['--', '-', '--'])

結果は以下のようになった。

このように、曲線でしか引けないような境界線も引くことができた。また特筆すべきは、楕円等のきれいな円ではなく、いびつな形で境界線を引けるということである。これにより、様々な形の分類に利用することが可能になる。