【Kaggle】タイタニックの分析をやってみる

Kaggleに登録してみたので備忘録を。
初心者ということもあり、まずは王道のタイタニックデータを使った分析をしてみた。

データの前処理


データはKaggle内にあり、TrainとTestにあらかじめ分けられている。
生存有無を予測する2値分類問題であり有名なデータ。

import numpy as np
import pandas as pd
from scipy import stats
 
#-- import
train = pd.read_csv("../input/titanic/train.csv")
test = pd.read_csv("../input/titanic/test.csv")


このデータは欠損があるようなので、今回は諸先輩方のブログを参考に中央値と最頻値で補完する。 また、カテゴリ変数は数値化しておく。

###--- 欠損処理
train["Age"] = train["Age"].fillna(train["Age"].median()) #--中央値で補完
train["Embarked"] = train["Embarked"].fillna("S") #-- 最頻値で補完 train["Embarked"].mode()

test["Age"] = test["Age"].fillna(test["Age"].median()) #--中央値で補完
test["Fare"] = test["Fare"].fillna(test["Fare"].median()) #--中央値で補完
 
###--- カテゴリの数値化
train["Sex"] = train["Sex"].map({'male':0, 'female':1})
train["Embarked"] = train["Embarked"].map({'S':0, 'C':1, 'Q':2})

test["Sex"] = test["Sex"].map({'male':0, 'female':1})
test["Embarked"] = test["Embarked"].map({'S':0, 'C':1, 'Q':2})


※あとから気づいたが、Embarked はダミー変数化した方が良かった。
pd.get_dummies(train, columns=['Embarked'])

バリデーションデータの準備


Kaggleでは、テストデータ(正解ラベルなし)があらかじめ準備されているため、次の手順で分析してみる。
* 訓練データを訓練&検証データに分割
* F1-スコアでモデル評価
* ベストモデルを用いてテストデータで予測

from sklearn import model_selection

Xval = ["Pclass", "Sex", "Age", "Fare", "Embarked"]

df_y = train["Survived"].values
df_X = train[Xval].values

train_X, valid_X, train_y, valid_y = model_selection.train_test_split(df_X, df_y, random_state=19)

test_X = test[Xval].values


モデル作成


初めてなので色々なモデルを試してみた(グリッドサーチ、ランダムサーチなどはしていない)。
順番に決定木、ランダムフォレスト、ロジスティック回帰(L2正則化、L1正則化、両方0.5ずつ)、線形サポートベクターマシーン(L2正則化、L1正則化)、サポートベクターマシーン(カーネル:Linear、rbf、poly、sigmoid)。

from sklearn.tree import DecisionTreeClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.svm import LinearSVC, SVC
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import classification_report, confusion_matrix, f1_score
 
modfnc = [DecisionTreeClassifier(random_state=19)
          ,RandomForestClassifier(random_state=19)
          ,LogisticRegression(random_state=19, penalty="l2")
          ,LogisticRegression(random_state=19, penalty="l1", solver='liblinear')
          ,LogisticRegression(random_state=19, penalty="elasticnet", solver='saga', l1_ratio=0.5)
          ,LinearSVC(random_state=19)
          ,LinearSVC(random_state=19, penalty='l1', loss='squared_hinge', dual=False)
          ,SVC(kernel='linear', random_state=19)
          ,SVC(kernel='rbf', random_state=19)
          ,SVC(kernel='poly', random_state=19)
          ,SVC(kernel='sigmoid', random_state=19)
         ]
 
score=0
 
for modfnc1 in modfnc:
    model = modfnc1
    model.fit(train_X, train_y)
    pred_y = model.predict(valid_X)
    print("----",modfnc1,"-----")
    #print(model.score(valid_X, valid_y))
    print(f1_score(valid_y, pred_y, average="micro"))
    #print(confusion_matrix(valid_y, pred_y))
    #print(classification_report(valid_y, pred_y))
    #print()
    
    if score < f1_score(valid_y, pred_y, average="micro"):
        best_model = modfnc1
        score = f1_score(valid_y, pred_y, average="micro")
 
print("---- best model ----")
print(best_model)
print(score)


このうち、F1-スコアでの評価ではランダムフォレストが最も良好だった。

---- DecisionTreeClassifier(random_state=19) -----
0.7713004484304933
---- RandomForestClassifier(random_state=19) -----
0.8430493273542601
---- LogisticRegression(random_state=19) -----
0.8071748878923767
---- LogisticRegression(penalty='l1', random_state=19, solver='liblinear') -----
0.8251121076233184
---- LogisticRegression(l1_ratio=0.5, penalty='elasticnet', random_state=19,
                   solver='saga') -----
0.6905829596412556
---- LinearSVC(random_state=19) -----
0.7309417040358744
---- LinearSVC(dual=False, penalty='l1', random_state=19) -----
0.8251121076233184
---- SVC(kernel='linear', random_state=19) -----
0.8116591928251121
---- SVC(random_state=19) -----
0.6636771300448431
---- SVC(kernel='poly', random_state=19) -----
0.6278026905829597
---- SVC(kernel='sigmoid', random_state=19) -----
0.5919282511210763
 
---- best model ----
RandomForestClassifier(random_state=19)
0.8430493273542601


テストデータで予測


ベストモデル(ランダムフォレスト)を使ってテストデータで予測を行う。
Kaggleのテストデータには正解ラベルが入っておらず、予測結果をKaggleの指定ページに提出してスコアを出す仕組みである(カンニングできない)。

best_model.fit(train_X, train_y)
pred_y2 = best_model.predict(test_X)
 
submission = pd.DataFrame({'PassengerId':test['PassengerId'], 'Survived':pred_y2})
submission.to_csv('submission.csv', index = False)


実際に結果を提出してみたところ、Score : 0.75119(44,476位)だった。

ハイパーパラメータを何もいじらずだったので、次はグリッドサーチやランダムサーチを使ってベストモデルを探してみたい。
もしくはポピュラーな LightGBMXGboost の実装をしてみたい(こっちの興味の方が強い)。

参考


【Kaggle超初心者向け】Titanicにチャレンジしてみた - Qiita
Kaggleコンペの提出データをKernelから提出する - Qiita
機械学習初心者がKaggleのTitanic課題でモデルを作る | 4番は司令塔

当ブログの関連記事:
【統計】ロジスティック回帰分析 - こちにぃるの日記
【機械学習_Python】Ridge回帰、Lasso回帰、Elastic Net回帰 - こちにぃるの日記
【機械学習_Python】決定木とランダムフォレスト - こちにぃるの日記
【機械学習_Python】サポートベクターマシーン - こちにぃるの日記

【SAS】call execute

executeの使い方の備忘録。

やりたいこと:
動的に生成したコードを実行したい。

call execute(コード) でいける。

data _null_;
  call execute("data a; a=1; run;");
run;


マクロ変数を投入することもできる。

%let CODE1 = '
    proc sql;
        create table a (
        col1 char 20 label="Columns 1"
        ,col2 num 8 label="Columns 2" format=8.1
        );
    quit;
';
data _null_;
  call execute(&CODE1.);
run;


execute内でコードを組み立てることもできる。

%let CODE2 = '
             col1 char 20 label="Columns 1"
            ,col2 num 8 label="Columns 2" format=8.1
    ';
data _null_;
  call execute(cats('proc sql; create table a (', &CODE., '); quit;'));
run;

【SAS】retainステートメント

文字列に関する記事はあまりなさそうだったので備忘録がてら。
個人メモで解説しないです。

やりたいこと:
* 文字型 ⇒ 複数行を1行にまとめる。
* 数値型 ⇒ 累積値を計算する。

サンプルデータ

data a;
  input a $ b;
  cards;
AA 1
BB 2
CC 3
;
run;


コード

data b;
  set a end=end;
  length a2 $200;
  retain a2 "" b2 0;      /* 初期値 */
  a2 = catx(" ", a2, a);  /* 文字列:空白スペースで連結 */
  b2 = b2 + b;            /* 数値型:累積 */
  if end then do;         /* 最終行をマクロ変数にする */
    call symputx("txt1", a2);
    call symputx("txt2", b2);
  end;
run;


output

%put "&txt1.";
 "AA BB CC"
 
%put "&txt2.";
 "6"

【機械学習_Python】k-means法

k-means法(k-平均法)についてメモ。

【目次】


前回取り上げた「k-近傍法」とは異なるので注意。
【機械学習_Python】k-近傍法(k-NN) - こちにぃるの日記

k-means法


データのクラスタリングに用いられる手法である。
分類問題における「教師なし学習」に属する。

アルゴリズムは比較的単純で、次の手順となる。

  1. 特徴空間上にk個のランダムな代表点(セントロイド)をプロットする。クラスタ数のkはハイパーパラメータであり、あらかじめ指定する。
  2. 各データとセントロイドの距離を求め、各データを最も近いセントロイドのグループに分類する。
  3. 上記2のグループでの重心(平均)を求め、新たなセントロイドとする。
  4. 上記2、3を繰り返す。
  5. セントロイドが更新されなくなれば終了し、分類を確定する。


f:id:cochineal19:20210524003450p:plain

数式としては、以下の関数を最小化する問題である。
「各クラスタ  C_{k} に属するデータ  x_{i} とセントロイド  \overline{x_{k}} との距離(ユーグリッド距離)の総和」(クラスタ内誤差平方和)を最小化できるセントロイドを求める。

クラスタ内誤差平方和(Sum of Squared errors of prediction、SSE)

 \quad J=\sum ^{K}_{k=1}\sum _{i\in C_{k}}\left\| x_{i}-\overline{x_{k}}\right\| ^{2}

 \qquad K:\verb|クラスタ数|,\ C_{k}:\verb|k番目のクラスタ|

 \qquad x_{i}:\verb|i番目のデータ|,\ \overline{x_{k}}:\verb|k番目のクラスタのセントロイド|


k-means++法


k-means法では、初期値のセントロイドをランダムな位置にプロットするため、複数のセントロイドが近くに位置する可能性があり、分類がうまくいかないことがある。

この問題の改良としてk-means++法があり、セントロイド同士が離れた位置にプロットされやすくなるアルゴリズムである。

エルボー法によるクラスタ数の推定


k-means法でのクラスタ数はハイパーパラメータであり、あらかじめ指定する必要がある。 このクラスタ数の推定方法として「エルボー法」がある。

エルボー法は、クラスタ数ごとのクラスタ内誤差平方和(SSE)を横に並べ、プロットがゆるやかになる位置を最適なクラスタ数とする方法である。

f:id:cochineal19:20210524021551p:plain:w450

Pythonコード


sklearn の cluster の KMeans で実装できる。
パラメタ説明はコメントアウトに記載する。

#--  k-means法
km1 = KMeans(n_clusters=4      # クラスター数指定。
             ,init="k-means++" # 'k-means++' or 'random'(k-means) 。default:'k-means++'
             ,n_init=10        # k-meansの実行回数
             ,random_state=0   # 乱数シード
             )
Y_km = km.fit_predict(X)  # 各データのクラスタ番号を求める
 
#-- SSE値を出力
print("Distortion: %.2f" % km.inertia_)


Irisデータでエルボー法のプロットを作ってみる。
クラスタ数(n_clusters=i)を1~10までループして、各クラスタ数のSSE(km.inertia_ で取得可)をプロットする。

from sklearn import datasets, cluster
import numpy as np
import matplotlib.pyplot as plt
 
#-- データ取得
X = datasets.load_iris()
 
#-- エルボー法(クラスタ数10まで)
SSE = []
for i in range(1, 11):
    km = cluster.KMeans(n_clusters = i
                ,init="k-means++"
                ,n_init=10
                ,random_state=1)
    km.fit(X.data)
    SSE.append(km.inertia_)
 
# グラフのプロット
plt.plot(range(1, 11), SSE, marker="o")
plt.xticks(np.arange(1, 11, 1))
plt.xlabel("クラスタ数")
plt.ylabel("クラスタ内誤差平方和(SSE)")
plt.show()

f:id:cochineal19:20210524025442p:plain:w300

参考


sklearn.cluster.KMeans — scikit-learn 0.24.2 documentation https://www.mext.go.jp/content/20200609-mxt_jogai01-000007843_007.pdf
https://www.jstage.jst.go.jp/article/fss/27/0/27_0_55/_pdf

【機械学習_Python】k-近傍法(k-NN)

k-近傍法(k-nearest neighbor algorithm、k-NN)についてメモ。
同じような名前の「k-means法」とは異なるので注意。

【目次】

k-近傍法とは


特徴空間において、距離が近い既知データから未知データのクラスを予測する方法である。 分類問題の「教師あり学習」である。
k-近傍法は教師データから学習(モデル作成)する訳でなく、単にデータ間の距離で分類する単純なアルゴリズムであり「怠惰学習」とも呼ばれる。

例えば、校庭(特徴空間)の様々な遊具(座標)の中で、すべり台付近に集まっている子たちが10人いて、 そのうち7人はA組、2人はB組の生徒だと分かっているとする(既知の情報)。 この情報から残り1名(未知の情報)を予測すると、A組の生徒の可能性が高いと予測される。

k-近傍法は、このようなアルゴリズムで周囲の既知データ(訓練データ)を多数決して、未知データのクラスを予測する。
なお、k-近傍法の「k」は周囲のいくつの既知データで予測するかを指定するハイパーパラメータである。多数決で予測を出すため奇数にするのが一般的である。
パラメータ k の値の違いで予測結果が変わることがある。

f:id:cochineal19:20210523124358p:plain:w450

なお、k=1 の場合を「最近傍法」と呼び、一番近くの既知データのクラスをそのまま採用する。

ミンコフスキー距離


ここで、近さの尺度として、特徴空間上の距離を表す「ミンコフスキー距離(Lpノルム)」に触れたい。
ミンコフスキー距離は「マンハッタン距離(L1ノルム)」や「ユーグリッド距離(L2ノルム)」を一般化したもので次の形をとる。

ミンコフスキー距離(Lpノルム)

\quad d\left( x,y\right) =\left( \sum ^{n}_{i=1}\left| x_{i}-y_{i}\right| ^{p}\right) ^{1/p}


パラメータ p=1 にてマンハッタン距離(L1ノルム)

\quad d\left( x,y\right) =\left( \sum ^{n}_{i=1}\left| x_{i}-y_{i}\right| ^{1}\right) ^{1/1} = \sum ^{n}_{i=1}\left| x_{i}-y_{i}\right|


パラメータ p=2 にてユーグリッド距離(L2ノルム)

\quad d\left( x,y\right) =\left( \sum ^{n}_{i=1}\left| x_{i}-y_{i}\right| ^{2}\right) ^{1/2} = \sqrt{ \sum ^{n}_{i=1}\left( x_{i}-y_{i}\right)^{2}}


k-近傍法では、このパラメータ p を設定することでマンハッタン距離やユーグリッド距離など複数の距離を扱える。一般的にはユーグリッド距離を使うことが多い。

pythonコード


sklearn の neighbors の KNeighbors にてk-近傍法を実装できる。
予測に用いる周囲の既知データ数(ハイパーパラメータ k )は、n_neighbors で指定する。デフォルトは n_neighbors=5
ミンコフスキー距離(Lpノルム)は p で指定する。デフォルトは p=2(ユーグリッド距離、L2ノルム)。

from sklearn import model_selection, neighbors, metrics
import seaborn as sns
import matplotlib.pyplot as plt

#-- データ取得・分割
iris = sns.load_dataset("iris")
X_train, X_test, y_train, y_test = train_test_split(iris.iloc[:,:4], iris.iloc[:,4], test_size=0.2, random_state=99)

#-- k-近傍法
kn = neighbors.KNeighborsClassifier(n_neighbors=5)
kn.fit(X_train, y_train)    #-- 学習
y_pred1 = kn.predict(X_test) #-- 予測

#-- 評価
print(metrics.confusion_matrix(y_test, y_pred1))
print(metrics.classification_report(y_test, y_pred1))


参考


sklearn.neighbors.KNeighborsClassifier — scikit-learn 0.24.2 documentation
https://www.mext.go.jp/content/20200609-mxt_jogai01-000007843_007.pdf
ミンコフスキー距離 - Wikipedia

【機械学習_Python】主成分分析(PCA)

主成分分析(Principal Component Analysis、PCA)についてのメモ

【目次】

主成分分析とは


相関のある特徴量から、互いに無相関の合成変数を作る手法である。この合成変数を主成分(Principal components)という。
主成分は特徴量の数だけ作ることができ、p次元の特徴量からp個の主成分ができる。このうち寄与率の高い主成分を選択することで、情報量を保ったまま次元削減することができる。例えば13次元の特徴量から2個の主成分を取り出し、2次元データとする。

f:id:cochineal19:20210521134924p:plain:w300

なお、主成分分析では特徴量の分散共分散行列を使う方法と相関行列を使う方法がある。元データが基準化(平均0、分散1)されていれば結果は同じになる。
特徴量間で単位が異なったり、分散の大きさが異なる場合は相関行列を使う方法が良い。

主成分軸の作り方


簡単のため、特徴量  x_{1}, x_{2} の2次元データ(特徴量行列X)で考える。
主成分分析では、特徴量 X を重み係数 W の方向に射影した主成分 Z の軸を作成する。

\quad \verb|行列表現: | Z=W^{T}X
\quad \verb|第1主成分: | z_{1}=w_{1}x_{1} + w_{2}x_{2}

次の概念図の通り、射影した主成分軸上での分散が大きいほど元の特徴量の情報を保つことができるため、主成分軸の作成は『主成分 Z の分散を最大化する重み係数 W を見つける問題』となる。

f:id:cochineal19:20210522231219p:plain:w500

この際、重み係数 W に制約がないと \infty になってしまうため、慣習的に重み係数 W のノルム(長さ)を1にして計算する(重み係数は割合が分かればよい)。

 \quad ||W|| = \sqrt{w_{p1}^{2} + w_{p2}^{2} + ... + w_{pp}^{2}} = 1


分散の最大化(ラグランジュの未定乗数法)

第1主成分を作る場合

特徴量  x_{1}, x_{2} の2次元データ(特徴量行列X)から第1主成分を作ることを考える。
つまり、 \sqrt{w_{1}^{2} + w_{2}^{2}} = 1 の制約の下、第1主成分  z_{1} の分散  \sigma_{z_{1}}^{2} の最大化を考える。


先ず、第1主成分の分散は、元データの分散と共分散に重み係数を掛けた形に書き換えられる。

 \quad \sigma_{z_{1}}^{2} = \dfrac{1}{n} \sum \left( z_{1} - \bar{z_{1}} \right)^{2}

 \qquad \ = \dfrac{1}{n} \sum \{ \left(w_{1}x_{1} + w_{2}x_{2}\right) - \left(w_{1}\overline{x_{1}} + w_{2}\overline{x_{2}} \right) \}^{2}

 \qquad \ = \dfrac{1}{n} \sum \{ w_{1}\left(x_{1} + \overline{x_{1}}\right) - w_{2}\left(x_{2} + \overline{x_{2}} \right) \}^{2}

 \qquad \ = w_{1}^{2}\dfrac{1}{n}\sum \left( x_{1}-\overline{x_{1}}\right)^{2} + w_{2}^{2}\dfrac{1}{n}\sum \left( x_{2}-\overline{x}_{2}\right) ^{2} + 2w_{1}w_{2} \dfrac{1}{n}\sum \left( x_{1}-\overline{x}_{1}\right) \left( x_{2}-\overline{x}_{2}\right)

 \qquad \ = w_{1}^{2} \sigma _{x_{1}}^{2} + w_{2}^{2} \sigma _{ x_{2}}^{2} + 2w_{1}w_{2} \sigma _{  x_{1} x_{2} }


これを行列表現で一般化すると、重みベクトル W と分散共分散行列 V から次の通り記載でき、重み係数 W を求める問題は W^{T}VW を最大化する問題へと帰着する。

 \quad \sigma_{z_{1}}^{2} = W^{T}VW = \begin{bmatrix} w_{1} & \cdot \cdot \cdot & w_{p}\end{bmatrix} \begin{bmatrix} \sigma _{  x_{1} x_{1} } & \cdot \cdot \cdot & \sigma _{  x_{p} x_{1} } \\ \vdots & \ddots & \vdots \\ \sigma _{  x_{1} x_{p} } & \ldots & \sigma _{  x_{p} x_{p} } \end{bmatrix}\begin{bmatrix} w_{1} \\ \vdots \\ w_{p} \end{bmatrix}


ここに制約条件 ||W||=1 を設け、ラグランジュ関数(  \lambdaラグランジュ乗数)を導入し、

 \quad L\left( w_{1},w_{2},\lambda \right) = \left( w_{1}^{2} \sigma _{x_{1}}^{2} + w_{2}^{2} \sigma _{ x_{2}}^{2} + 2w_{1}w_{2} \sigma _{  x_{1} x_{2} } \right) - \lambda\left( w_{1}^{2}+w_{2}^{2}-1 \right)

 \qquad \qquad \qquad \quad \ =W^{T}VW- \lambda\left( W^{T}W-1 \right)

 w_{1}, w_{2}, \lambda偏微分より次を得る。

 \quad \dfrac{\partial L}{w_{1}}=2w_{1}\sigma _{x_{1}}^{2}+2w_{2}\sigma _{x_{1}x_{2}}-2w_{1}\lambda=w_{1}\sigma _{x_{1}}^{2}+w_{2}\sigma _{x_{1}x_{2}}-w_{1}\lambda=0

 \quad \dfrac{\partial L}{w_{2}}=2w_{2}\sigma _{x_{2}}^{2}+2w_{1}\sigma _{x_{1}x_{2}}-2w_{2}\lambda=w_{2}\sigma _{x_{1}x_{2}}+w_{1}\sigma _{x_{2}}^{2}-w_{2}\lambda=0

 \quad \dfrac{\partial L}{\partial \lambda }=-\left( w_{1}^{2}+w_{2}^{2}-1\right) =0


この重み係数  w_{1}, w_{2} の解について、変形して並べてみる。

 \quad w_{1}\sigma _{x_{1}}^{2}+w_{2}\sigma _{x_{1}x_{2}}=w_{1} \lambda

 \quad w_{1}\sigma _{x_{1}x_{2}}+w_{2}\sigma _{x_{2}}^{2}=w_{2} \lambda

これを行列表現で一般化すると次のようになる。

 \quad VW=W \lambda I

 \quad \begin{bmatrix} \sigma _{x_{1}x_{1}} & \ldots & \sigma _{x_{p} x_{1}} \\  \vdots & \ddots & \vdots \\ \sigma _{x_{1} x_{p}} & \ldots & \sigma _{x_{2}x_{2}} \end{bmatrix}\begin{bmatrix} w_{1} \\  \vdots \\ w_{p} \end{bmatrix}=\begin{bmatrix} w_{1} \\  \vdots \\ w_{p} \end{bmatrix}\begin{bmatrix} \lambda &  \ldots & 0  \\  \vdots & \ddots & \vdots \\ 0 & \ldots  & \lambda \end{bmatrix}


最後の式は、分散共分散行列  V固有値問題であり、固有値固有ベクトルはそれぞれ次に対応する。
・最大固有値  \lambda =主成分  z_{1} の分散  \sigma_{z_{1}}^{2}
固有ベクトル=重み係数  W

つまり、重み係数 W を求める問題は、分散共分散行列  V固有値問題に帰着する。
なお、標準化する場合、分散共分散行列でなく相関行列 R を用いる。


第p主成分を作る場合

第2主成分以上を作りたい場合も基本的な計算は変わらず、最終的には分散共分散行列  V固有値問題に帰着する。

 \quad VW=W \lambda I

 \quad \begin{bmatrix} \sigma _{x_{1}x_{1}} & \ldots & \sigma _{x_{p} x_{1}} \\  \vdots & \ddots & \vdots \\ \sigma _{x_{1} x_{p}} & \ldots & \sigma _{x_{2}x_{2}} \end{bmatrix}\begin{bmatrix} w_{1} \\  \vdots \\ w_{p} \end{bmatrix}=\begin{bmatrix} w_{1} \\  \vdots \\ w_{p} \end{bmatrix}\begin{bmatrix} \lambda &  \ldots & 0  \\  \vdots & \ddots & \vdots \\ 0 & \ldots  & \lambda \end{bmatrix}


なお、第2主成分は第1主成分、第3主成分は第1・第2主成分と直交する軸のうち分散を最大化する軸を求める。以降も同様。

ここで、固有値を大きい順に  \lambda_{1}, \lambda_{2}, ..., \lambda_{p} と並べ、対応する固有ベクトル w_{}, w_{2}, ... , w_{p} とすると、この並び順が第1主成分、第2主成分、・・第p主成分となる。先述の通り、最大固有値  \lambda _{1} が第1主成分である。


次元削減の基準

寄与率と累積寄与率

第k主成分でどのくらいの分散を表せているか(情報量があるか)を固有値を用いて「寄与率」という尺度で表すことができる。

 \quad \dfrac{\lambda_{k}}{\sum _{i=1}^{p} \lambda} = \dfrac{\lambda_{k}}{\lambda_{1}+\lambda_{2}+...+\lambda_{p}}

また、第1主成分から第k主成分までの寄与率の合計を「累積寄与率」という。

 \quad \dfrac{\sum _{j=1}^{k} \lambda_{j}}{\sum _{i=1}^{p} \lambda_{i}} = \dfrac{\lambda_{1}+\lambda_{2}+...+\lambda_{k}}{\lambda_{1}+\lambda_{2}+...+\lambda_{k}+...+\lambda_{p}}


この「累積寄与率」を十分に高く保てる範囲で次元削減する方法がある。一般的に80%、50%などあるが根拠があるわけではない。

カイザー基準とスクリー基準

分散の加法性により、基準化したデータセット(平均0、分散1)における分散の合計は、特徴量の数と等しくなる(2次元なら分散の合計は2、p次元なら分散の合計はp)。よって、固有値が1に満たないものは、特徴量1個分の情報量もないと考えることができる。このように固有値が1以上であるかで判断する方法を「カイザー基準」という。

また、固有値を大きい順に並べてプロットし、推移がゆるやかになるところまで選択する方法がある。この方法を「スクリー基準」と言い、この時の図をスクリープロットと言う。

f:id:cochineal19:20210523020003p:plain:w400

計算手順のまとめ


主成分分析~次元削減までの手順は、簡単に次の通りまとめられる。
①:p次元の特徴量行列 X を標準化(平均0、分散1)する。
②:①から相関行列 R を求める。
③:②から固有値  \lambda 固有ベクトル(重み係数  W )を求める。
④:③の固有値  \lambda を大きい順にk個取り出し、これらに対応する固有ベクトル(重み係数  W )を結合して特徴量の変換行列 D を作る。
⑤:④の変換行列 D を用いて①の特徴量行列 X を変換する。X'=XD

Pythonコード


sklearn の decomposition にある PCA で主成分分析ができる。
次元削減したい場合は PCA の n_components で指定する。例:2次元データに圧縮する。PCA(n_components=2)

from sklearn.decomposition import PCA
import numpy as np

#-- 標準化
X = (X - X.mean(axis=0)) / X.std(axis=0)

#-- 主成分分析
pca1   = PCA() #-- 次元数を指定する場合、n_components
X_pca1 = pca1.fit_transform(X) #-- 分析

print("-- 固有値(主成分の分散) --")
print(pca1.explained_variance_)
print("-- 固有ベクトル(重み係数W) --")
print(pca1.components_)
print("-- 寄与率 --")
print(pca1.explained_variance_ratio_)


参考


統計学(出版:東京図書), 日本統計学会編
高橋 信, 井上 いろは, マンガでわかる統計学 [因子分析編], トレンドプロ
市川 伸一 他, SASによるデータ解析入門 [第3版], 東京大学出版会
sklearn.decomposition.PCA — scikit-learn 0.24.2 documentation
http://www.ie.reitaku-u.ac.jp/~tak/datB/datB_prin00.pdf
http://manabukano.brilliant-future.net/document/text-PCA.pdf
http://web.wakayama-u.ac.jp/~wuhy/am10.pdf

【機械学習_Python】サポートベクターマシーン

SVM(サポートベクターマシーン)についてのメモ。


SVM(サポートベクターマシーン)とは


SVMは「マージン最大化」の考えに基づいて分類する手法である。
n次元空間に超平面(線形識別関数、二次元なら直線)で分類境界を作ることでベクトルを分離する。
(マージンが大きいほど、未知データでの識別誤りが起こりにくいだろうという理屈でマージン最大化する)。

手法としては、線形分離可能を仮定する「ハードマージン法」と任意の識別誤りを許容する「ソフトマージン法」がある。
さらに、非線形問題に対しても、高次元特徴空間を作り「カーネルトリック」を使って分類する方法がある。
前者を線形SVM、後者を非線形SVMと呼ぶ。

線形SVM(ハードマージン法とソフトマージン法)


「ハードマージン法」は、訓練データ(学習するデータ)を超平面で完全に分離できると仮定した問題である。

f:id:cochineal19:20210520010026p:plain

ただし、現実的には線形分離可能な問題は珍しく、ある程度の識別誤りを許容する「ソフトマージン法」が用いられる。識別誤りのパラメータとしてスラック変数 ξ を導入し、ハイパーパラメータCでその許容を調整する。

f:id:cochineal19:20210520010101p:plain

主問題と双対問題

SVMでは「主問題」を直接解くのではなく、計算コストの良い「双対問題」に帰着して解く方法が実装されている。双対問題を解く際、後に記載する「カーネルトリック」も扱える。

双対問題を解くにあたっては、主問題を不等式制約のラグランジュ未定乗数法を用いて双対問題へ置き換える(計算の詳細は割愛)。
以下はハードマージン法とソフトマージン法の式を横に並べたもの。赤字部分が異なる箇所。
f:id:cochineal19:20210520021551p:plain

ここで双対問題とは、主変数  w, b, \xi を最小化し、双対変数  \alpha , \beta を最大化する問題である。
これは主変数にとっての最小値が、双対変数にとっての最大値になる「鞍点」(あんてん)を探している(鞍点が最適解ということ。鞍点はKKT条件を満たす)。鞍点は字のごとく、馬の鞍(くら)の形。

※双対問題・鞍点の概念は、以下リンクのスライド22枚目が分かりやすかった(他サイト様のもの)。
はじめてのパターン認識 第8章 サポートベクトルマシン

以上により最適解 \alpha_{i} が求まれば、ラグランジュ関数にて求めた次の式で係数wが求まる。

 \quad w=\sum_{i\in sv}α_{i}t_{i}x_{i}, \  sv:サポートベクトル

ここで、マージンの外側(本記事の概念図で背景色を付けた領域)にある正解ベクトルは \alpha_{i}=0 となり、マージンの内側(  0 < \xi _{i} \leq 1, または  1 < \xi_{i} )にあるベクトルは  \alpha _{i} = C となる。そして、ちょうどマージンの上(概念図の点線)にあるベクトルは  0 < \alpha_{i} < C となることから、\alpha ≠ 0 であるマージンの内側とマージン上のベクトル(これらが「サポートベクトル」と言われる)のみを用いて係数wを求めれば良い。


<再掲>
f:id:cochineal19:20210520010101p:plain:w400

Pythonによる実行

sklearn の svm.LinearSVC で線形SVMが実行できる。
ソフトマージン法のハイパーパラメータCはデフォルトで C=1.0 である。

Pythonコード

from sklearn.svm import LinearSVC
 
mod1 = LinearSVC(C=C, random_state=1)
mod1.fit(X_train, y_train)
y_pred1 = mod1.predict(X_test)
print(confusion_matrix(y_test, y_pred1))
print(classification_report(y_test, y_pred1))
print()


なお、線形SVMでもL1正則化とL2正則化を指定できる。デフォルトは penalty="l2"(L2正則化)。
正則化については次の記事を参照。
【機械学習_Python】Ridge回帰、Lasso回帰、Elastic Net回帰 - こちにぃるの日記

非線形SVMカーネルトリック


ソフトマージン法によって線形分離不可能なデータセットも扱えるようになったが、現実的には非線形の分類境界を作りたい場面が多い。
このような場面に対して、入力空間を高次元特徴空間に変換(特徴量を非線形写像)して分類可能とする方法がある。
f:id:cochineal19:20210520162542p:plain:w400

上の例では、関数  \phi \left( x_{i}\right)^{T} \phi \left( x_{j}\right) を用いて新たな特徴量(例えば  x_{1} x_{2} )を作り、3次元空間とすることで分離可能になった。このように、特徴量を高次元に写像するほど分離可能性が高まる。
この方法は、次元を高めるほど計算コストが爆発的に高まり実装しづらい欠点があったが、カーネル関数を導入した「カーネルトリック」と呼ばれる方法によって計算コスト面を解消することができる。

カーネル関数

 \quad K\left(x_{i},x_{j}\right)=\phi \left( x_{i}\right)^{T} \phi \left( x_{j}\right)

※簡単な記載になってしまうが、カーネル関数を用いることで関数  \phi \left( x_{i}\right)^{T} \phi \left( x_{j}\right) の中身を計算せずとも、等価な結果を計算することができる。

f:id:cochineal19:20210520163237p:plain

Pythonによる実行

sklearn の svm.SVC非線形SVMが実行できる。
ソフトマージン法と同じくハイパーパラメータ:Cを設定可能。デフォルトは C=1.0
カーネルの指定は、linearpolyrbfsigmoidprecomputed があり、デフォルトは kernel='rbf'poly多項式カーネルrbf がGaussカーネルsigmoid がシグモイドカーネルで、linear は線形SVMsklearn.svm.LinearSVC と少し異なるらしい)。
多クラス分類においては、decision_function_shapeにおいて、ovo(one vs one、1対1)とovr(one vs rest、1対残り)を指定できる。ovoはクラス同士のペアを総当たりするため計算量が多い。デフォルトは decision_function_shape='ovr'

Pythonコード

from sklearn.svm import SVC
 
mod1 = SVC(C=C, kernel='rbf', random_state=1)
mod1.fit(X_train, y_train)
y_pred1 = mod1.predict(X_test)
print(confusion_matrix(y_test, y_pred1))
print(classification_report(y_test, y_pred1))
print()


Wineデータで実行


Wineデータから線形分離が難しそうなデータの組を適当に見つけて、線形SVM非線形SVMの識別性能を比較してみる。

import seaborn as sns
import pandas as pd
import numpy as np
from sklearn import datasets, model_selection
from sklearn.svm import LinearSVC, SVC
from sklearn.metrics import classification_report, confusion_matrix
import matplotlib.pyplot as plt

ROW = datasets.load_wine()
ADS = pd.DataFrame(ROW.data, columns=ROW.feature_names).assign(target=np.array(ROW.target))
sns.pairplot(ADS[["ash","proline","target"]], hue="target")

f:id:cochineal19:20210520192333p:plain:w400

特徴量 "ash", "proline" の target 1 vs. 2 が良い感じに混ざっていたので、今回これらを用いる。

ADS12 = ADS[ADS["target"]!=0]
ADS12["target"] = ADS["target"] - 1
 
X_train, X_test, y_train, y_test = model_selection.train_test_split(
        ADS12[["ash","proline"]], ADS12["target"], test_size=0.3, random_state=999)
 
#-- 線形SVM
mod1 = LinearSVC(C=1, random_state=999)
mod1.fit(X_train, y_train)
y_pred1 = mod1.predict(X_test)
print('線形SVM')
print(confusion_matrix(y_test, y_pred1))
print(mod1.score(X_test, y_test))
#print(classification_report(y_test, y_pred1))
print()
 
#-- 非線形SVM
klist=['poly','rbf','sigmoid']
for k in klist:
    mod1 = SVC(C=1, kernel=k, random_state=999)
    mod1.fit(X_train, y_train)
    y_pred1 = mod1.predict(X_test)
    print('非線形SVM ', k)
    print(confusion_matrix(y_test, y_pred1))
    print(mod1.score(X_test, y_test))
    #print(classification_report(y_test, y_pred1))
    print()

f:id:cochineal19:20210520192618p:plain:w150

線形SVMと比べ、テストデータに対する accuracy は、非線形SVM多項式カーネル(poly)とGaussカーネル(rbf)が良さそうだった。

参考


https://home.hiroshima-u.ac.jp/tkurita/lecture/svm.pdf
サポートベクターマシン(Support Vector Machine, SVM)~優秀な(非線形)判別関数~
サポートベクターマシンの詳しい理論的な解説について【線形分離可能な場合】
はじめてのパターン認識 第8章 サポートベクトルマシン

本ブログは個人メモです。 本ブログの内容によって生じた損害等の一切の責任を負いかねますのでご了承ください。