機械学習レポート

目次

• 機械学習
  • 第一章　線形回帰モデル
  • 第二章　非線形回帰モデル
  • 第三章　ロジスティック回帰モデル
  • 第四章　主成分分析
  • 第五章　アルゴリズム
  • 第七章　サポートベクターマシン

別視点からのレポート

統計学レポート
データ分析と知識発見レポート

機械学習とは

• タスクTを性能指標Pで測定し、その性能が経験Eにより改善される場合、タスクTおよび性能指標Pに関して経験Eから学習すると言われている

  • タスクT：アプリケーション
  • 経験E：データ

  （トム・ミッチェル　1997）

• 人がプログラムするのは認識の仕方ではなく学習の仕方（数学で記述）

線形・非線形

線形：直線
非線形：曲線（曲がり角がある）

入力・出力

入力：説明変数、または、特徴量（m次元ベクトル）
出力：目的変数（スカラー）

データ分割

未学習・過学習

性能評価

線形回帰モデル：平均二乗誤差（残差平方和）、最小二乗法
非線形回帰モデル：ホールドアウト法、クロスバリデーション（交差検証）、グリッドサーチ

第一章　線形回帰モデル

単回帰：説明変数が一次元（m=1）、直線
重回帰：説明変数が多次元（m>1）、曲線

###【実装演習結果】

単回帰

訓練データ作成

n_sample = 100
var = .2

def linear_func(x):
    return 2 * x + 5

def add_noise(y_true, var):
    return y_true + np.random.normal(scale=var, size=y_true.shape)

def plt_result(xs_train, ys_true, ys_train):
    plt.scatter(xs_train, ys_train, facecolor="none", edgecolor="b", s=50,  label="training data")
    plt.plot(xs_train, ys_true, label="$2 x + 5$")
    plt.legend()

#データの作成
xs = np.linspace(0, 1, n_sample)
ys_true = linear_func(xs)
ys = add_noise(ys_true, var)

print("xs: {}".format(xs.shape))
print("ys_true: {}".format(ys_true.shape))
print("ys: {}".format(ys.shape))

#結果の描画
plt_result(xs, ys_true, ys)

学習

１次関数 y(x)=ax+b における、 a と b を求める。
訓練データ X=[x1,x2,...,xn]T,y=[y1,y2,...,yn]T に対して、最小化する目的関数は L=∑ni=1(yi−(axi+b))2 と書け、
∂L∂a=−2∑ni=1(yi−(axi+b))xi=0
∂L∂b=−2∑ni=1(yi−(axi+b))=0⋯(1)
より、目的関数を最小にする a,b は以下のように求まる。
(∑ni=1x2i∑ni=1xi∑ni=1xin)(a^b^)(a^b^)===(∑ni=1xiyi∑ni=1yi)1n∑ni=1x2i−(∑ni=1xi)2(n−∑ni=1xi−∑ni=1xi∑ni=1x2i)(∑ni=1xiyi∑ni=1yi)1n∑ni=1x2i−(∑ni=1xi)2(n∑ni=1xiyi−∑ni=1xi∑ni=1yi∑ni=1x2i∑ni=1yi−∑ni=1xi∑ni=1xiyi)⋯(2)
(1), (2)から a^=Cov[x,y]/Var[x],b^=μy−a^μx で求める。

def train(xs, ys):
    cov = np.cov(xs, ys, ddof=0)
    a = cov[0, 1] / cov[0, 0]
    b = np.mean(ys) - a * np.mean(xs)
    return cov, a, b

cov, a, b = train(xs, ys)
print("cov: {}".format(cov))
print("coef: {}".format(a))
print("intercept: {}".format(b))

cov: [[0.08501684 0.16491662]
[0.16491662 0.35375915]]
coef: 1.9398113591765453
intercept: 5.070786959564312

予測

入力に対する値を y(x)=ax+b で予測する

xs_new = np.linspace(0, 1, n_sample)
ys_pred = a * xs_new + b

plt.scatter(xs, ys, facecolor="none", edgecolor="b", s=50, label="training data")
plt.plot(xs_new, ys_true, label="$2 x + 5$")
plt.plot(xs_new, ys_pred, label="prediction (a={:.2}, b={:.2})".format(a, b))
plt.legend()
plt.show()

多項式回帰

訓練データ生成

n_sample = 10
var = .25

def sin_func(x):
    return np.sin(2 * np.pi * x) 

def add_noise(y_true, var):
    return y_true + np.random.normal(scale=var, size=y_true.shape)

def plt_result(xs, ys_true, ys):
    plt.scatter(xs, ys,facecolor="none", edgecolor="b", s=50,  label="training data")
    plt.plot(xs, ys_true, label="$\sin(2\pi x)$")
    plt.legend()

#データの作成
xs = np.linspace(0, 1, n_sample)
ys_true = sin_func(xs)
ys = add_noise(ys_true, var)

print("xs: {}".format(xs.shape))
print("ys_true: {}".format(ys_true.shape))
print("ys: {}".format(ys.shape))

#結果の描画
plt_result(xs, ys_true, ys)

学習

モデルとして以下を用いる。
y(x)=∑di=0wixi=wTϕ(x)
ただし、 w=[w0,w1,...,wd]T,ϕ(x)=[1,x,x2,...,xd]T である。
訓練データ X,y に対しては y=Φw と書ける。
ただし、 Φ=[ϕ(x1),ϕ(x2),..,ϕ(xn)]T である。
よって、最小化する目的関数は L=||y−Φw||2 と書け、
∂L∂w=−2ΦT(y−Φw)=0 より、求める回帰係数 w は以下のように書ける。
w^=(ΦTΦ)−1ΦTy

def polynomial_features(xs, degree=3):
    """多項式特徴ベクトルに変換
         X = [[1, x1, x1^2, x1^3],
                 [1, x2, x2^2, x2^3],
                   ...                             
                 [1, xn, xn^2, xn^3]]"""
    X = np.ones((len(xs), degree+1))
    X_t = X.T #(100, 4)
    for i in range(1, degree+1):
        X_t[i] = X_t[i-1] * xs
    return X_t.T

Phi = polynomial_features(xs)
Phi_inv = np.dot(np.linalg.inv(np.dot(Phi.T, Phi)), Phi.T)
w = np.dot(Phi_inv, ys)

予測

入力を多項式特徴ベクトル ϕ(x) に変換し、 y=w^ϕ(x) (y(x)=Φw^) で予測する。

Phi_test = polynomial_features(xs)
ys_pred = np.dot(Phi_test, w)

plt.scatter(xs, ys, facecolor="none", edgecolor="b", s=50, label="training data")
plt.plot(xs, ys_true, label="$\sin(2\pi x)$")
plt.plot(xs, ys_pred, label="prediction")
# for i in range(0, 4):
#     plt.plot(xs, Phi[:, i], label="basis")
plt.legend()
plt.show()

重回帰分析

訓練データ生成 (3次元入力)

np.random.random((10, 3))

n_sample = 100
var = .2

def mul_linear_func(x):
    global ww # 追加
    ww = [1., 0.5, 2., 1.]
    return ww[0] + ww[1] * x[:, 0] + ww[2] * x[:, 1] + ww[3] * x[:, 2]

def add_noise(y_true, var):
    return y_true + np.random.normal(scale=var, size=y_true.shape)

def plt_result(xs_train, ys_true, ys_train):
    plt.scatter(xs_train, ys_train, facecolor="none", edgecolor="b", s=50,  label="training data")
    plt.plot(xs_train, ys_true, label="$2 x + 5$")
    plt.legend()

x_dim = 3

X = np.random.random((n_sample, x_dim))
ys_true = mul_linear_func(X)
ys = add_noise(ys_true, var)

学習

モデルとして以下を用いる。
y(x)=∑di=0wixi=wTx
ただし、陽には書かないが、 x には定数項のための 1 という要素があることを仮定する。
訓練データ X,y に対しては y=Xw と書ける
よって、最小化する目的関数は L=||y−Xw||2 と書け、
∂L∂w=−2XT(y−Xw)=0 より、求める回帰係数 w は以下のように書ける
w^=(XTX)−1XTy

def add_one(x):
    return np.concatenate([np.ones(len(x))[:, None], x], axis=1)

X_train = add_one(X)
# pinv = np.dot(np.linalg.inv(np.dot(X_train.T, X_train)), X_train.T)
# w = np.dot(pinv, y_train)

予測

入力に対する値を y(x)=w^Tx (y=Xw^) で予測する

パラメータ推定結果

for i in range(len(w)):
    print("w{0}_true: {1:>5.2}   w{0}_estimated: {2:>5.2}".format(i, ww[i], w[i]))

w0_true: 1.0 w0_estimated: -0.49
w1_true: 0.5 w1_estimated: 1.3e+01
w2_true: 2.0 w2_estimated: -3.4e+01
w3_true: 1.0 w3_estimated: 2.2e+01

###【考察】
###【演習問題】
　線形回帰は教師あり学習の回帰手法の1つであり、目的変数が説明変数の係数に対して線形なモデルを学習するものである。線形回帰では訓練データに適合する回帰係数を学習する際、二重誤差の最小化を行う。

第二章　非線形回帰モデル

• 基底展開法

  • 回帰関数として、基底関数と呼ばる既知の非線形関数とパラメータベクトルの線形結合を使用
  • 未知パラメータは線形回帰モデルと同様に最小2乗法や最尤法により推定

• よく使われる基底関数

  • 多項式関数
  • ガウス型基底関数
  • スプライン関数/Bスプライン関数

• 正則化法（罰則化法）

  • 「モデルの複雑さに伴って、その値が大きくなる正則化項（罰則項）を課した関数」を最小化
  • 正則化項（罰則項）
    • 形状によっていくつもの種類があり、そろぞれ推定量の性質が異なる
    • L1ノルム：Ridge 推定量
    • L2ノルム：Lasso 推定量

###【実装演習結果】

#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)

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)

#Ridge

from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Ridge

kx = rbf_kernel(X=data, Y=data, gamma=50)
#KX = rbf_kernel(X, x)

#clf = LinearRegression()
clf = Ridge(alpha=30)
clf.fit(kx, target)

p_ridge = clf.predict(kx)

plt.scatter(data, target,label='data')
for i in range(len(kx)):
    plt.plot(data, kx[i], color='black', linestyle='-', linewidth=1, markersize=3, label='rbf', alpha=0.2)

#plt.plot(data, p, color='green', marker='o', linestyle='-', linewidth=0.1, markersize=3)
plt.plot(data, p_ridge, color='green', linestyle='-', linewidth=1, markersize=3,label='ridge regression')
#plt.legend()

print(clf.score(kx, target))

from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline

#PolynomialFeatures(degree=1)

deg = [1,2,3,4,5,6,7,8,9,10]
for d in deg:
    regr = Pipeline([
        ('poly', PolynomialFeatures(degree=d)),
        ('linear', LinearRegression())
    ])
    regr.fit(data, target)
    # make predictions
    p_poly = regr.predict(data)
    # plot regression result
    plt.scatter(data, target, label='data')
    plt.plot(data, p_poly, label='polynomial of degree %d' % (d))

#Lasso

from sklearn.metrics.pairwise import rbf_kernel
from sklearn.linear_model import Lasso

kx = rbf_kernel(X=data, Y=data, gamma=5)
#KX = rbf_kernel(X, x)

#lasso_clf = LinearRegression()
lasso_clf = Lasso(alpha=10000, max_iter=1000)
lasso_clf.fit(kx, target)

p_lasso = lasso_clf.predict(kx)

plt.scatter(data, target)

#plt.plot(data, p, color='green', marker='o', linestyle='-', linewidth=0.1, markersize=3)
plt.plot(data, p_lasso, color='green', linestyle='-', linewidth=3, markersize=3)

print(lasso_clf.score(kx, target))

from sklearn import model_selection, preprocessing, linear_model, svm

# SVR-rbf
clf_svr = svm.SVR(kernel='rbf', C=1e3, gamma=0.1, epsilon=0.1)
clf_svr.fit(data, target)
y_rbf = clf_svr.fit(data, target).predict(data)
 
# plot

plt.scatter(data, target, color='darkorange', label='data')
plt.plot(data, y_rbf, color='red', label='Support Vector Regression (RBF)')
plt.legend()
plt.show()

from sklearn.model_selection import train_test_split
x_train, x_test, y_train, y_test = train_test_split(data, target, test_size=0.1, random_state=0)

from keras.callbacks import EarlyStopping, TensorBoard, ModelCheckpoint

cb_cp = ModelCheckpoint('/content/drive/My Drive/study_ai_ml/skl_ml/out/checkpoints/weights.{epoch:02d}-{val_loss:.2f}.hdf5', verbose=1, save_weights_only=True)
cb_tf  = TensorBoard(log_dir='/content/drive/My Drive/study_ai_ml/skl_ml/out/tensorBoard', histogram_freq=0)

def relu_reg_model():
    model = Sequential()
    model.add(Dense(10, input_dim=1, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='relu'))
    model.add(Dense(1000, activation='linear'))
#     model.add(Dense(100, activation='relu'))
#     model.add(Dense(100, activation='relu'))
#     model.add(Dense(100, activation='relu'))
#     model.add(Dense(100, activation='relu'))
    model.add(Dense(1))

    model.compile(loss='mean_squared_error', optimizer='adam')
    return model

from keras.models import Sequential
from keras.layers import Input, Dense, Dropout, BatchNormalization
from keras.wrappers.scikit_learn import KerasRegressor

# use data split and fit to run the model
estimator = KerasRegressor(build_fn=relu_reg_model, epochs=100, batch_size=5, verbose=1)

history = estimator.fit(x_train, y_train, callbacks=[cb_cp, cb_tf], validation_data=(x_test, y_test))

y_pred = estimator.predict(x_train)

plt.title('NonLiner Regressions via DL by ReLU')
plt.plot(data, target, 'o')
plt.plot(data, true_func(data), '.')
plt.plot(x_train, y_pred, "o", label='predicted: deep learning')
#plt.legend(loc=2)

print(lasso_clf.coef_)

[-0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
 -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0. -0.
 -0. -0. -0. -0. -0. -0. -0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.
  0.  0.  0.  0.  0.  0.  0.  0.  0.  0.]

###【考察】
###【演習問題】

第三章　ロジスティック回帰モデル

• 分類問題
  • 目的変数が 0 or 1 の値
• 入力
  • 実数
• 出力
  • シグモイド関数により、 0 or 1
• 分類の評価方法
  • 再現率（Recall）
  • 適合率（Precision）
• F値
• 誤差関数
  • 交差エントロピー（クロエントロピー）誤差関数

###【実装演習結果】

訓練データ生成

n_sample = 100
harf_n_sample = 50
var = .2

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

def plt_data(x_train, y_train):
    plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train, facecolor="none", edgecolor="b", s=50, label="training data")
    plt.legend()

#データ作成
x_train, y_train = gen_data(n_sample, harf_n_sample)
#データ表示
plt_data(x_train, y_train)

識別モデルとして p(y=1|x;w)=σ(wTx) を用いる。
ただし、 σ(⋅) はシグモイド関数であり、 σ(h)=11+exp(−h) で定義される。
また、陽には書かないが、 x には定数項のための 1 という要素があることを仮定する。

学習

訓練データ X=[x1,x2,...,xn]T,y=[y1,y2,...,yn]T(yi={0,1}) に対して尤度関数 L は以下のように書ける。
L(w)=∏ni=1p(yi=1|xi;w)yi(1−p(yi=1|xi;w))1−yi
負の対数尤度関数は
−logL(w)=−∑ni=1[yilogp(yi=1|xi;w)+(1−yi)log(1−p(yi=1|xi;w))]
のように書ける。 これを最小化する w を求める。
dσ(h)dh=σ(h)(1−σ(h)) と書けることを考慮し、負の対数尤度関数を w で偏微分すると、
∂∂w(−logL(w))==−∑i=1n[yi(1−σ(wTxi))−(1−yi)σ(wTxi)]xi∑i=1n(σ(wTxi)−yi))xi
この式が 0 となる w は解析的に求められないので、今回は −logL(w) の最小化問題を最急降下法を用いて解く。
最急降下法では学習率を η とすると、以下の式で w を更新する。
w←w−η∂∂w(−logL(w))

def add_one(x):
    return np.concatenate([np.ones(len(x))[:, None], x], axis=1)

def sigmoid(x):
    return 1 / (1 + np.exp(-x))

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

X_train = add_one(x_train)
max_iter=100
eta = 0.01
w = sgd(X_train, max_iter, eta)

予測

入力に対して、 y=1 である確率を出力する。よって
p(y=1|x;w)=σ(wTx) の値が
0.5 より大きければ1に、小さければ0に分類する。

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 = add_one(xx)
proba = sigmoid(np.dot(X_test, w))
y_pred = (proba > 0.5).astype(np.int)

plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(xx0, xx1, proba.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))

###【考察】
###【演習問題】
ロジスティック回帰は教師あり学習の分類手法の1つであり、2値分類では、入力が一方のクラスに属する確率を出力するように学習する。ロジスティック回帰では学習の際、尤度関数の最大化を行う。

• 確率的勾配降下法
  • 訓練データサイズが大きい場合はバッチ最急降下法より確率的勾配降下法を用いるのが良い
  • 確率的勾配降下法を適用するために目的関数のパラメータに関する勾配が計算できる必要がある
  • 訓練データの順番を変えると結果も変わる可能性がある

第四章　主成分分析（PCA、Principle Component Analysis）

• 寄与率：第k主成分の全分数に対する割合（第k主成分が持つ情報量の割合）
  • 第1~元次元分の主成分で表示した時、固有値の和と元のデータの分散が一致
• 累積寄与率：第1-k主成分まで圧縮した際の情報損失量の割合

###【実装演習結果】

訓練データ生成

n_sample = 100

def gen_data(n_sample):
    mean = [0, 0]
    cov = [[2, 0.7], [0.7, 1]]
    return np.random.multivariate_normal(mean, cov, n_sample)

def plt_data(X):
    plt.scatter(X[:, 0], X[:, 1])
    plt.xlim(-5, 5)
    plt.ylim(-5, 5)

X = gen_data(n_sample)
plt_data(X)

学習

訓練データ X=[x1,x2,...,xn]T に対して E[x]=0 となるように変換する。
すると、不偏共分散行列は Var[x]=1n−1XTX と書ける。
Var[x] を固有値分解し、固有値の大きい順に対応する固有ベクトルを第1主成分( w1 ), 第2主成分( w2 ), ...とよぶ。

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

def plt_result(X, first, second):
    plt.scatter(X[:, 0], X[:, 1])
    plt.xlim(-5, 5)
    plt.ylim(-5, 5)
    # 第１主成分
    plt.quiver(0, 0, first[0], first[1], width=0.01, scale=6, color='red')
    # 第２主成分
    plt.quiver(0, 0, second[0], second[1], width=0.01, scale=6, color='green')

#分散共分散行列を標準化
meean, stan_cov = get_moments(X)
#固有値と固有ベクトルを計算
eigenvalues, eigenvectors = np.linalg.eigh(stan_cov)
components = get_components(eigenvectors, n_components)

plt_result(X, eigenvectors[0, :], eigenvectors[1, :])

変換（射影）

元のデータを m 次元に変換(射影)するときは行列 W を W=[w1,w2,⋯,wm] とし、データ点 x を z=WTx によって変換(射影)する。
よって、データ X に対しては Z=XTW によって変換する。

def transform_by_pca(X, pca):
    mean = X.mean(axis=0)
    return np.dot(X-mean, components)

Z = transform_by_pca(X, components.T)
plt.scatter(Z[:, 0], Z[:, 1])
plt.xlim(-5, 5)
plt.ylim(-5, 5)

Z = transform_by_pca(X, components.T)
plt.scatter(Z[:, 0], Z[:, 1])
plt.xlim(-5, 5)
plt.ylim(-5, 5)

逆変換

射影されたデータ点 z を元のデータ空間へ逆変換するときは x¯=(WT)−1z=Wz によって変換する。
よって、射影されたデータ Z に対しては X¯=ZWT によって変換する。

mean = X.mean(axis=0)
X_ = np.dot(Z, components.T) + mean

plt.scatter(X_[:, 0], X_[:, 1])
plt.xlim(-5, 5)
plt.ylim(-5, 5)

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(X)

print('components: {}'.format(pca.components_))
print('mean: {}'.format(pca.mean_))
print('covariance: {}'.format(pca.get_covariance())) 

plt_result(X, pca.components_[0, :], pca.components_[1, :])

from sklearn.decomposition import PCA
pca = PCA(n_components=2)
pca.fit(X)
plt_result(X, pca.components_[0, :], pca.components_[1, :])

print('components: {}'.format(pca.components_))
print('mean: {}'.format(pca.mean_))
print('covariance: {}'.format(pca.get_covariance()))

###【考察】
　ディープラーニングのオートエンコーダーでも次元圧縮する。
###【演習問題】

主成分分析は教師なし学習の1つであり、データから重要な成分をみつける手法である。主成分におけて重要な成分とはデータの分散の大きい成分であり、この成分のことを主成分という。各主成分は互いに直交するように選ばれる。

• 主成分分析の使用目的は以下である

• データの特徴を抽出する

• 高次元データを低次元にする

• 高次元データを可視化する

• 主成分の数は以下である。

• 主成分の数の最大値はデータの特徴量の数である

• 一般に可視化の際は、主成分の数は3以下にすることが多い

• 主成分の数が少ないほど、データを射影した際に情報が落ちる

第五章　アルゴリズム

• アルゴリズム
• 問題を解くといった作業を行う場合の処理の手順のこと

k近傍法（kNN、k-Nearest Neighbor）

• 教師なし学習
• クラス分類
• 分類のための機械学習手法
  • 最近傍のデータをk個取ってきて、それらがもっとも多く所属するクラスに識別
  • kを変化させると結果も変わる
  • k=1 のときは、最近傍法という

###【実装演習結果】

訓練データ生成

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)

予測するデータ点との、距離が最も近い 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 = 3

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)

k平均法（k-means、k-means clustering）

• 教師なし学習
• クラスタリング手法
• 与えられたデータをk個のクラスタに分類する
• k平均法（k-means）のアルゴリズム
  1. 各クラスタ中心の初期値を設定する
  2. 各データに対して、各クラスタ中心との距離を計算し、最も距離が近いクラスタを割り当てる
  3. 各クラスタの平均ベクトル（中心）を計算する
  4. 収束するまで 2, 3の処理を繰り返す
• 中心の初期値を変えるとクラスタリング結果も変わりうる

###【実装演習結果】

訓練データ生成

def gen_data():
    x1 = np.random.normal(size=(100, 2)) + np.array([-5, -5])
    x2 = np.random.normal(size=(100, 2)) + np.array([5, -5])
    x3 = np.random.normal(size=(100, 2)) + np.array([0, 5])
    return np.vstack((x1, x2, x3))

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

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

def plt_result(X_train, centers, xx):
    #　データを可視化
    plt.scatter(X_train[:, 0], X_train[:, 1], c=y_pred, cmap='spring')
    # 中心を可視化
    plt.scatter(centers[:, 0], centers[:, 1], s=200, marker='X', lw=2, c='black', edgecolor="white")
    # 領域の可視化
    pred = np.empty(len(xx), dtype=int)
    for i, x in enumerate(xx):
        d = distance(x, centers)
        pred[i] = np.argmin(d)
    plt.contourf(xx0, xx1, pred.reshape(100, 100), alpha=0.2, cmap='spring')

y_pred = np.empty(len(X_train), dtype=int)
for i, x in enumerate(X_train):
    d = distance(x, centers)
    y_pred[i] = np.argmin(d)

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

plt_result(X_train, centers, xx)

N予備校 機械学習 入門 機械学習 k-means

k平均法++（k-means++）

• k平均法に加えて初期値をランダムにするように変える

N予備校 機械学習 入門 機械学習 k-means++

###【考察】
###【演習問題】

• k-means は教師なしクラスタリング手法の1つである。
• k-means のアルゴリズムにおいて最適な解を得るために用いられる工夫
• クラスタ中心の初期値を変えて実行する
• 繰り返し回数を調整する
• クラスターの数を調整する
• k-means と k-NNはいずれのアルゴリズムも距離計算を
  行う

第七章　サポートベクターマシン

サポートベクターマシン(SVM, support vector machine)

• 2クラス分類のための機械学習手法
• 線形モデルの正負で2値分類

###【実装演習結果】

訓練データ生成① （線形分離可能）

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)

学習

特徴空間上で線形なモデル y(x)=wϕ(x)+b を用い、その正負によって2値分類を行うことを考える。
サポートベクターマシンではマージンの最大化を行うが、それは結局以下の最適化問題を解くことと同じである。
ただし、訓練データを X=[x1,x2,...,xn]T,t=[t1,t2,...,tn]T(ti={−1,+1}) とする。
minw,b12||w||2
subject toti(wϕ(xi)+b)≥1(i=1,2,⋯,n)
ラグランジュ乗数法を使うと、上の最適化問題はラグランジュ乗数 a(≥0) を用いて、以下の目的関数を最小化する問題となる。
L(w,b,a)=12||w||2−∑ni=1aiti(wϕ(xi)+b−1)⋯(1)
目的関数が最小となるのは、 w,b に関して偏微分した値が 0 となるときなので、
∂L∂w=w−∑ni=1aitiϕ(xi)=0
∂L∂b=∑ni=1aiti=aTt=0
これを式(1)　に代入することで、最適化問題は結局以下の目的関数の最大化となる。
L~(a)==∑i=1nai−12∑i=1n∑j=1naiajtitjϕ(xi)Tϕ(xj)aT1−12aTHa
ただし、行列 H の i 行 j 列成分は Hij=titjϕ(xi)Tϕ(xj)=titjk(xi,xj) である。また制約条件は、 aTt=0(12||aTt||2=0) である。
この最適化問題を最急降下法で解く。目的関数と制約条件を a で微分すると、
dL~da=1−Ha
dda(12||aTt||2)=(aTt)t
なので、 a を以下の二式で更新する。
a←a+η1(1−Ha)
a←a−η2(aTt)t

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)

予測

新しいデータ点 x に対しては、 y(x)=wϕ(x)+b=∑ni=1aitik(x,xi)+b の正負によって分類する。
ここで、最適化の結果得られた ai(i=1,2,...,n) の中で ai=0 に対応するデータ点は予測に影響を与えないので、 ai>0 に対応するデータ点（サポートベクトル）のみ保持しておく。 b はサポートベクトルのインデックスの集合を S とすると、 b=1S∑s∈S(ts−∑ni=1aitik(xi,xs)) によって求める。

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=['--', '-', '--'])

# マージンと決定境界を可視化
plt.quiver(0, 0, 0.1, 0.35, width=0.01, scale=1, color='pink')

訓練データ生成② （線形分離不可能）

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)

学習

元のデータ空間では線形分離は出来ないが、特徴空間上で線形分離することを考える。

今回はカーネルとして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=['--', '-', '--'])

ソフトマージンSVM

訓練データ生成③（重なりあり）

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)

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

学習

分離不可能な場合は学習できないが、データ点がマージン内部に入ることや誤分類を許容することでその問題を回避する。

スラック変数 ξi≥0 を導入し、マージン内部に入ったり誤分類された点に対しては、 ξi=|1−tiy(xi)| とし、これらを許容する代わりに対して、ペナルティを与えるように、最適化問題を以下のように修正する。

minw,b12||w||2+C∑ni=1ξi

subject toti(wϕ(xi)+b)≥1−ξi(i=1,2,⋯,n)

ただし、パラメータ C はマージンの大きさと誤差の許容度のトレードオフを決めるパラメータである。この最適化問題をラグランジュ乗数法などを用いると、結局最大化する目的関数はハードマージンSVMと同じとなる。

L~(a)=∑ni=1ai−12∑ni=1∑nj=1aiajtitjϕ(xi)Tϕ(xj)

ただし、制約条件が ai≥0 の代わりに 0≤ai≤C(i=1,2,...,n) となる。（ハードマージンSVMと同じ ∑ni=1aiti=0 も制約条件）

X_train = x_train
t = np.where(y_train == 1.0, 1.0, -1.0)

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

C = 1
eta1 = 0.01
eta2 = 0.001
n_iter = 1000

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.clip(a, 0, C)

予測

index = a > 1e-8
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(-4, 4, 100), np.linspace(-4, 4, 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=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=['--', '-', '--'])

N予備校 機械学習 入門 機械学習 サポートベクターマシン

###【考察】
###【演習問題】