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

線形回帰モデル

・教師あり学習、回帰
・回帰問題とは、m次元ベクトルの入力から、スカラ値（目的変数）を予測する
・線形は比例（直線、平面、超平面）。変数１つは単回帰、複数だと重回帰
$$\hat{y}=\sum_{j=1}^m w_j x_j + w_0 + 誤差$$　（wがパラメータ）

＊データの分割（学習用データと検証用データ）による汎化性能測定に留意

・最小二乗法
　平均二乗誤差（MSE：真値と予測値の二乗誤差の平均）が最小となるものを探索
　＊外れ値に弱いためデータの前処理が重要
　MSEを最小とするためMSEをwに関して微分したものが0となる点を求める
　$$MSE=\frac{1}{n}\sum_{i=1}^n (\hat{y_i}-y_i)^2$$
$$\frac{\partial}{\partial w}MSE=0$$

・【実習】sklearnとボストンデータセット（skl_regression.ipynb）
　目的：ボストンの住宅価格の予測

＊データセットの中を見て、上限切り捨てのデータや外れ値等の不自然なデータがないか確認

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

# 説明変数（部屋数）
data = df.loc[:, ['RM']].values
# 目的変数
target = df.loc[:, 'PRICE'].values

# 予測
model.predict([[1]])

array([-25.5685118])

＊データを部屋数だけで線形回帰（単回帰）すると、１部屋の物件がマイナス価格に。学習データは６部屋ぐらいのものなので、１部屋の予想は難しい（外挿問題には弱い）

# 説明変数（犯罪率、部屋数）の重回帰
data2 = df.loc[:, ['CRIM', 'RM']].values

# 予測
model2.predict([[0.2, 6]])

array([21.04870738])

model2.predict([[0.8, 6]])
array([20.88975943])

犯罪率を変えてもさほど価格は変わらない。部屋数が大きく影響

非線形回帰モデル

・教師あり学習、回帰
・基底関数として非線形関数を使用：$x_i$だったところが$f(x_i)$になる
・よく使われる基底関数
　多項式関数：$x^j$
　ガウス型基底関数：$exp{\frac{(x-\mu_j)^2}{2h_j}}$
・線形回帰モデルと同様に最小二乗法による
・未学習（underfitting）：モデルの表現力が低くデータを表現できない
　過学習（overfitting）：学習データに最適化しすぎ、汎化性能がでない
　　→学習データを増やす
　　　不要な基底関数の削除
　　　正規化法：誤差関数に複雑さに対する罰則項を加える
　　　　　　　（ex.$MSE+\lambda ww^T$）
正規化パラメータはハイパーパラメータ
・検証方法（学習データと検証データの分割
　ホールドアウト法：学習と検証のためにデータを分け、変更しない。データが少ないと精度が上がらない。
　クロスバリデーション（交差検証）：ｋ個に分割し、１つを検証用とし評価。検証用を変えｋ回評価して、平均の制度をとる。
・ハイパーパラメータの選択
　グリッドサーチ：すべてのハイパーパラメータの組み合わせを評価し決定
　ベイズ最適化（BayesianOptimization）：すべてのパラメータを試すのはコストが高いため、良さそうなところを選んで確認する。

・【実習】sklearnとボストンデータセット（skl_nonlinear regression.ipynb）
　真の関数（$z = 1-48x+218x^2-315x^3+145x^4$）にノイズを加えデータを生成

# カーネルリッジ回帰
from sklearn.kernel_ridge import KernelRidge
# Radial basis function kernelを使用。ハイパーパラメータ0.0002
clf = KernelRidge(alpha=0.0002, kernel='rbf')

# Radial basis function kernelを使用。ハイパーパラメータ0.1
clf = KernelRidge(alpha=0.1, kernel='rbf')

https://scikit-learn.org/stable/modules/generated/sklearn.kernel_ridge.KernelRidge.html
ハイパーパラメータαは正規化の強さ。大きくすると緩やかな曲線となり過学習をより抑制する

ロジスティクス回帰モデル

・教師あり学習、２クラス分類（回帰とついているが、分類）、識別的アプローチ
・分類問題とは、m次元ベクトルの入力から、0か1（目的変数）を出力する
・アプローチの種類
　識別的アプローチ：$p(C_k | x)$を直接モデル化
　生成的アプローチ：$p(C_k)$と$p(x|C_k)$をモデル化し、ベイズの定理により$p(C_k | x)$を求める。

・シグモイド関数
　０から１の間に変換し、確率として利用
　$$\sigma(z) = \frac{1}{1+e^(-z)}$$
$$\frac{\partial \sigma(z)}{\partial z}=\sigma(z)(1-\sigma(z))$$

・最尤推定
　データから最もらしいパラメータを推定
　尤度関数を最大化するパラメータ$w$を推定
$$P(y_1,y_2,…,y_n|w_0,w_1,…,w_m)$$

　＊計算を簡単にし、桁落ちを防ぐため尤度関数の対数をとる
　＊マイナスをつけて最小化するパラメータを求める。

・勾配降下法
　解析的に解くのが困難なため、傾きを求め下がる方向に最小値を探索
　学習率で傾きに応じてどの程度移動するかを調整　
$$w^(k+1)=w^k-\mu\sum_{i=1}^n (y_i-p_i)x_i$$

・確率的勾配降下法（ＳＧＤ）
　勾配降下法はn個の和を求める必要があるため、nが増えると必要メモリが増加
　ランダムにデータを選択してパラメータを更新

・実習（np_logistic_regression.ipynb）
　ビデオでは実習がカットされてる（？）
　実装を見るため、numpy版のコード中身を確認

# SGDの実装（データをひとつづつ学習し、wを更新）
def sgd(X_train, max_iter, eta):
    w = np.zeros(X_train.shape[1])
    for _ in range(max_iter):
        w_prev = np.copy(w)
        sigma = sigmoid(np.dot(X_train, w))
        grad = np.dot(X_train.T, (sigma - y_train))
        w -= eta * grad
        if np.allclose(w, w_prev):
            return w
    return w

max_iter：学習の反復の最大回数
eta：学習率

$p(y=1|x;w)=σ(w^Tx)$を確率とみて0.5より大きければ1、小さければ0

proba = sigmoid(np.dot(X_test, w))
y_pred = (proba > 0.5).astype(np.int)

・実習（skl_logistic_regression.ipynb）
　目的：タイタニックの乗客データからその生死を予測
　googleドライブからcsvを読み込むが、パス名が違っているので修正

　今回は、年齢と船室等級の２つの入力で学習

from sklearn.linear_model import LogisticRegression
model=LogisticRegression()
model.fit(data, label)

/usr/local/lib/python3.7/dist-packages/sklearn/utils/validation.py:760: DataConversionWarning: A column-vector y was passed when a 1d array was expected. Please change the shape of y to (n_samples, ), for example using ravel().
  y = column_or_1d(y, warn=True)
LogisticRegression(C=1.0, class_weight=None, dual=False, fit_intercept=True,
                   intercept_scaling=1, l1_ratio=None, max_iter=100,
                   multi_class='auto', n_jobs=None, penalty='l2',
                   random_state=None, solver='lbfgs', tol=0.0001, verbose=0,
                   warm_start=False)

ライブラリではデフォルトで正則化がされているのを確認。デフォルトの動作を要確認
https://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html

　年齢が高く、船室等級が低いほど死亡が多い傾向

主成分分析（ＰＣＡ）

　動画が４分ぐらいで終わっていろいろ不安になる。

・教師なし学習、データの次元削減
・次元削減により3次元以下することで可視化できたり、計算コストを下げることができる
　　可視化ができるが、そこからそれぞれの軸が何を意味するかを読み取るのは困難
・削減後の分散が最大になる（情報が保存される）ように第１主成分を決定し、さらに第１主成分と直行で分散が最大になるように次の主成分を決定していく
　→第１主成分を基準に新しい座標系に置き換える

・実習（skl_pca.ipynb）

　がんの検診データ。ロジスティクス回帰モデルでデータから悪性か良性かを分類している。
　ロジスティクス回帰モデルの項目でのタイタニックのデータの実習の際は、はじめから２つの変数を選択して学習したが、今回は３０個の変数をもとに学習しているため、このままでは表等にプロットして可視化することができない。

# PCA
# 次元数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]

　寄与率は、元のデータの散らばり具合をどれだけ残しているかを表す。
　今回は２次元に圧縮しているので、それぞれの次元の寄与率がでている。
　合わせて62なので、元のデータから４割弱の情報が落ちていると考えられる。
　一般的にはにはこの寄与率を見てどこまで圧縮するかを判断することになる。

　２次元に圧縮することにより、プロットして可視化することができる。

・実習（np_pca.ipynb）

　numpyでの実装を確認する。

n_components=2

def get_moments(X):
    mean = X.mean(axis=0)
    stan_cov = np.dot((X - mean).T, X - mean) / (len(X) - 1)
    return mean, stan_cov
    
def get_components(eigenvectors, n_components):
#     W = eigenvectors[:, -n_components:]
#     return W.T[::-1]
    W = eigenvectors[:, ::-1][:, :n_components]
    return W.T

　データから共分散行列を固有値展開し、固有値が大きいものから主成分としていく。
　ここでは２つを選択

アルゴリズム

　レポートの単元名が「アルゴリズム」なのは携帯電話を携帯と略すぐらい意味を失っている。
　講義のシラバス上では、「最近傍・K-近傍アルゴリズム」「K-meansアルゴリズム」が記載されているので、これらのことと考えられる。
　一般的には、「K-近傍法（K-NN）」「k平均法（K-means）」と記載するケースが多いと思われる。
　表記ポリシーがわからない。いろいろ不安になる。
　この２つは似ているため混同しがちであるが、全く別の概念。分類とクラスタリングの違いに留意

・K-近傍法（K-NN）
　教師あり学習、分類
　未知のデータについて、既知のデータから最も近いものを順にｋ個選択し、その多数決でクラスを決定する方法
　ｋ＝１のとき特に「最近傍法」と呼ばれる
　一般にｋが大きくなるほど境界面が滑らかになり、外れ値に強くなる

・K平均法（K-means）
　教師なし学習、クラスタリング
　K個のクラスタにクラスタリングする
　アルゴリズム
　　K個のクラスタの重心を求め、重心に近いものをそのクラスタに再分類することを繰り返すことでクラスタリングする
　クラスタの重心の初期値は最初にランダムで決定するが、初期値によって適切にクラスタリングされない場合がある。対策として、何度か初期値を変えて評価を繰り返す方法や、K-means++という初期重心が適度に離れて選定されるようにしたアルゴリズムがある

・実習（np_knn.ipynb）
　K-NNの実装を確認

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

　事前に訓練するような手順はなく、未知のデータを学習データと比較する

・実習（np_kmeans.ipynb）
　K-meansの実装を確認

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

n_clusters = 3
iter_max = 100

# 各クラスタ中心をランダムに初期化
centers = X_train[np.random.choice(len(X_train), n_clusters, replace=False)]

for _ in range(iter_max):
    prev_centers = np.copy(centers)
    D = np.zeros((len(X_train), n_clusters))
    # 各データ点に対して、各クラスタ中心との距離を計算
    for i, x in enumerate(X_train):
        D[i] = distance(x, centers)
    # 各データ点に、最も距離が近いクラスタを割り当
    cluster_index = np.argmin(D, axis=1)
    # 各クラスタの中心を計算
    for k in range(n_clusters):
        index_k = cluster_index == k
        centers[k] = np.mean(X_train[index_k], axis=0)
    # 収束判定
    if np.allclose(prev_centers, centers):
        break

　次元数の増加、訓練データの増加で計算コストが増えていくので留意する

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

　ついにビデオがなくなり参考図書の紹介になっている。
　サンプルコードがあるので確認する。

・教師あり学習、２クラス分類
・マージン（境界面から最も近い各クラスのデータ（これをサポートベクターという。）からの距離）が最大になるように境界面を決定する。
・線形分離可能な場合は、ハードマージンと呼ばれる。
・線形分離不可能な場合は、ソフトマージンと呼ばれる誤判定を許容して境界面を決定する。この際、誤差をペナルティとして加える（正則化）
・$y(x)=wϕ(x)+b$の正負で２値分類を行う

・線形分離できないデータについて、カーネル関数により高次元へ投影してから分離する方法をカーネル法という。
・一般的に、高次元への投影は計算コストが高いため、直接的に座標変換せずに内積の評価を行うことで計算量を削減する手法をカーネルトリックという。このカーネルトリックはＳＶＭのほかに主成分分析等さまざまな手法と組み合わせて利用できる。
　最終的に必要なのは新しい座標系での内積の和（評価値）であるので、これをもとの座標の関数と置き換えることで計算を削減する
・代表的なカーネル関数はRBFカーネル（ガウシアンカーネル）
　$K(x,x')=exp(-\gamma||x-x'||^2)$

・実習（np_svm.ipynb）

　numpyでの実装をみつつ、論理を確認
　境界の関数 $y(x)=wϕ(x)+b$とサポートベクトルの距離の最大化を考える

　目的関数はラグランジュ関数により$L(w,b,a)=\frac{1}{2}||w||2−∑_{i=1}^n a_it_i(wϕ(x_i)+b−1)$ （ラグランジュ乗数 a）。これを最小化する。

　a を以下の二式で更新する。

　$a←a+η\1(1−Ha) $

　$a←a−η_2(aT^t)t$

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)

eta1、eta2はハイパーパラメータ、n_iterは実施回数として与える。

　これによって得られた境界線の正負で分類する
　$y(x)=wϕ(x)+b=∑_{i=1}^n a_it_ik(x,x_i)+b$

　分類された様子を可視化したもの。

　次に、下記のような線形分離不可能なデータの場合で実施する

# RBFカーネル
def rbf(u, v):
        sigma = 0.8
        return np.exp(-0.5 * ((u - v)**2).sum() / sigma**2)

RBFカーネルにより座標系を高次元に投影し、線形分離可能な形にしたうえで学習させる。

　

　分類された様子を可視化したもの。

　最後に、外れ値等で完全に分類できない場合についてソフトマージンにより分類する。

　

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

　clipメソッドによりaが0からCの範囲に制約される。Cは正則化係数のハイパーパラメータでどれだけ誤分類にペナルティを加えるかを意味する。ハードマージンのコードと比較すると、Cが無限に大きい場合は、ハードマージンの時と同じ制約条件になることがわかる。