lazy diary

数学とかコンピュータサイエンスとか

非線形SVMによる分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

ハードマージンSVM、ソフトマージンSVMと続いて今回でSVMは最後です。
今までは全て線形SVMでしたが、今回は非線形SVMとなります。非線形SVMでは特徴空間次元の内積計算が多数必要になり、それを避けるためにカーネル法という手法が用いられます。カーネル法にも種類は様々あるのですが今回は最も一般的なガウスカーネルを用います。データはソフトマージンSVMのときと同じデータです。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

class NonLinearSoftMarginSVM():
  
  def __init__(self, X, y, lam, sigma):
    self.X = X
    self.y = y
    self.lam = lam
    self.sigma = sigma
    
  #ガウスカーネル
  def kernel(self, x, y):
    return np.exp(np.dot(x - y, x - y) / (-2 * self.sigma))
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * self.kernel(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    #条件緩和によって変化した部分
    G_1 = (-1) * np.identity(num_of_data)
    G_2 = np.identity(num_of_data)
    G = matrix(np.vstack((G_1, G_2)))
    h_1 = np.zeros((num_of_data, 1))
    h_2 = np.ones((num_of_data, 1)) * self.lam
    h = matrix(np.vstack((h_1, h_2)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    #2次計画問題終了
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon or alpha[i] > self.lam: #条件緩和によって変化した部分
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    print(Sup_index)
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      for j in range(num_of_data):
        b = b + alpha[j] * label.iloc[j] * self.kernel(self.X.iloc[j], self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index, alpha]

前回のソフトマージンSVMと比べて変化しているのは、内積計算が全てカーネルによる計算に置き換えられているという点だけです。続いて、ここにデータを入れてみます。

lam = 1.0
sigma = 1.0
X = data[["x", "y"]]
y = data["c"]
label = data["c"] * 2 - 1
svm = NonLinearSoftMarginSVM(X , y, lam, sigma)
[w, b, Sup_index, alpha] = svm.fit()

ここまでは実はさほど計算時間はかかりません。決定境界をプロットしてみましょう。

def f(x, alpha, Sup_index, b, X, label):
  result = b
  for i in range(len(Sup_index)):
    j = Sup_index[i]
    result = result + alpha[j] * label[j]  * svm.kernel(X.iloc[j], x)
  return result

#色の設定
col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")

#メッシュグリッドを作成して決定境界をプロット
seq = np.linspace(0, 10, 50)
X1, X2 = np.meshgrid(seq, seq)
Z = []
x1 = X1.reshape(-1, 1)
x2 = X2.reshape(-1, 1)
for i in range(X1.size):
  Z.append(f(np.array([x1[i][0], x2[i][0]]), alpha, Sup_index, b, X, label))

plt.title("Non Lnear Soft Margin SVM")
plt.xlabel("x")
plt.ylabel("y")
Z = np.array(Z).reshape(X1.shape)
plt.contour(X1, X2, Z, colors = "black", levels = [0.0])
plt.scatter(data["x"], data["y"], c = col)

f:id:mt19:20181022155837p:plain
このように非線形の決定境界となります。SVM自体の話はここまでで、あとはどうやってモデルを良いものにしていくかが難しいところです。例えば今回はソフトマージンなので正則化係数、さらにはガウスカーネルなので分散\sigmaを事前に設定する必要があります。他には多項式カーネルであれば2つのハイパーパラメータを設定しなければなりません。クロスバリデーションなどを適用しながらこれらを適切に設定することが求められます。

参考

非線形の決定境界を書く際に参考になりました。
python で等高線を描くなら meshgrid して contour | コード7区

ソフトマージンSVMによる分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

前回はハードマージンSVMを実装しました。今回はそれの拡張のようなもので、ソフトマージンSVMを行ないます。ハードマージンSVMではデータは完全に線形分離できるものにしてくださいと書きました。
ハードマージンSVMによる分類 - lazy diary
しかし現実問題そんな都合の良いデータばかりではないでしょう。なんらかの理由でラベルの異なるグループに入ってしまうデータもあるはずです。ソフトマージンSVMでは、そのようなデータが入ることを許すよう条件を緩和します。しかし、どの程度まで許すかということが問題になります。そこで条件を緩和しすぎず、なおかつマージンを大きくすることを目指したのがソフトマージンSVMです。
実際に解析するにあたってその程度を決める指標として正則化係数を私たちは与える必要があります。これについてはクロスバリデーションで決める方法がおそらく最もシンプルだと思われます。

コード自体はハードマージンSVMとほとんど変わりません。二次計画問題で、条件が増える点と、サポートベクトル集合にも条件がさらに加わることだけ注意します。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

class SoftMarginSVM():
  
  def __init__(self, X, y, lam):
    self.X = X
    self.y = y
    self.lam = lam
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * np.dot(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    #条件緩和によって変化した部分
    G_1 = (-1) * np.identity(num_of_data)
    G_2 = np.identity(num_of_data)
    G = matrix(np.vstack((G_1, G_2)))
    h_1 = np.zeros((num_of_data, 1))
    h_2 = np.ones((num_of_data, 1)) * self.lam
    h = matrix(np.vstack((h_1, h_2)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    #2次計画問題終了
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon or alpha[i] > self.lam: #条件緩和によって変化した部分
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    print(Sup_index)
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      b = b + np.dot(w, self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index]

データを作ります。ハードマージンSVMと違い多少は入り乱れたデータでも平気です。

#create data
x_1 = 4.2
y_1 = 7
x_2 = 7
y_2 = 4.2
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

これをSVMにかけてみましょう。

lam = 10.0
svm = SoftMarginSVM(data[["x", "y"]], data["c"], lam)
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.title("Soft Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022142912p:plain
このように、ハードマージンSVMの厳しすぎる条件を緩和することでソフトマージンSVMでは分類を達成します。

ハードマージンSVMによる分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

今回はハードマージンSVMというものを実装します。
ソフトでもない上に非線形でもないという最も基本的なSVMです。ハードマージンSVMさえ理解してしまえば他のSVMはその拡張みたいなものです。

必要なものをimportしておきます。

from matplotlib import pyplot as plt
import numpy as np
import pandas as pd
import cvxopt
from cvxopt import matrix

このcvxoptというのは二次計画問題のソルバーです。詳しくはこの記事の最後にあるブログを参考にしてください。google colabにもデフォルトではインストールされていないようなので、

pip install cvxopt

で導入しておきましょう。

データは相変わらずいい加減に作ったものです。とは言いつつも今回はハードマージンSVMですので、線形分離できないデータは困ります。かならず超平面で分離できるデータをチョイスしてください。2次元なら直線で完全に分離できるものです。

#create data
x_1 = 3
y_1 = 7
x_2 = 7
y_2 = 2
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

ここからが本編です。流れとしては
1. データから二次計画問題を解く
2. パラメータwを求める
3. サポートベクトル集合を求める
4. サポートベクトル集合から切片bを求める
となっています。

class HardMarginSVM():
  
  def __init__(self, X, y):
    self.X = X
    self.y = y
    
  def fit(self):
    num_of_data = self.X.shape[0]
    label = self.y * 2.0 - 1.0 #ラベルを-1と1に変換する
    
    #2次計画問題を解く
    alpha = np.zeros((num_of_data, 1))
    q = matrix(np.ones((num_of_data ,1)))
    P = np.zeros((num_of_data, num_of_data))
    for i in range(num_of_data):
      for j in range(num_of_data):
        P[i,j] = label.iloc[i] * label.iloc[j] * np.dot(self.X.iloc[i], self.X.iloc[j])
    P = matrix(P)
    A = matrix(np.array(label), (1, num_of_data))
    b = matrix(0.0, (1,1))
    G = matrix((-1) * np.identity(num_of_data))
    h = matrix(np.zeros((num_of_data, 1)))
    sol = cvxopt.solvers.qp(P, -q, G, h, A, b)
    alpha = sol["x"]
    
    #パラメータを求める
    w = 0
    for i in range(num_of_data):
      w = w + alpha[i] * label.iloc[i] * self.X.iloc[i]
    
    #切片bを求めるためのサポートベクター集合S
    epsilon = 0.00001 #これより小さいものを0とみなす
    Sup_index = []
    for i in range(num_of_data):
      if alpha[i] < epsilon:
        continue
      Sup_index.append(i)
    b = 0
    temp = 0
    for i in range(len(Sup_index)): #分離超平面より上と下の個数を数える
      if label.iloc[Sup_index[i]] == 1:
        temp = temp + 1
      b = b + np.dot(w, self.X.iloc[Sup_index[i]])
    b = ((temp - (len(Sup_index) - temp)) - b) / len(Sup_index)
    
    return [w, b, Sup_index]

これを上のデータに当てはめてみましょう。

svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022105408p:plain
黄色になっている部分はサポートベクトルです。SVMの特徴として実質的にサポートベクトルにしか依存しないので、次のような外れ値が入っても影響がありません。

data.loc[100] = [0, 20, 1]
svm = HardMarginSVM(data[["x", "y"]], data["c"])
[w, b, Sup_index] = svm.fit()

col = []
for i in range(len(data["c"])):
  if data["c"].iloc[i] == 1:
    col.append("red")
  else:
    col.append("blue")
seq = np.arange(np.min(data["x"]) - 1, np.max(data["x"]), 0.1)
plt.scatter(data["x"], data["y"], c = col)
plt.plot(seq, (-w[0] * seq - b) / (w[1]), c = "black")
plt.scatter(data["x"].iloc[Sup_index], data["y"].iloc[Sup_index], c = "yellow", marker = "D")
plt.title("Hard Margin SVM")
plt.xlabel("x")
plt.ylabel("y")

f:id:mt19:20181022110158p:plain

参考

今回、二次計画問題を解くにあたって参考にしました。SMOアルゴリズムについてはまだ理解できていないため、今回はcvxoptというソルバーに頼りました。
二次計画法(Quadratic Programming)の概要とPythonソルバーcvxoptの使い方 - MyEnigma

正準判別によるマルチクラス多次元データの可視化

参考書は相変わらずこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

今回は正準判別というものを行ないます。
目的は、可視化の難しい多次元データをなんとかして私たちの理解できる2 or 3次元に落とし込むことです。正準判別の考え方はフィッシャーの線形判別と似ていて、各グループに属するデータをグループ内の分散を小さく、かつ、グループ間の平均差を大きく取れる方向へ射影しようという考え方です。なるべくグループが分かれるように低次元に射影するということです。
これに似ているものとして主成分分析がありますが、これは多次元データの情報をなるべく崩さないように低次元に持っていく方法であり正準判別とは少々目的も方法も異なります。

さて、正準判別を実装していきます。多次元データを自分で作るのはあれなのでirisのデータを用います。これを可視化することが目的です。

from sklearn import datasets
import pandas as pd
import numpy as np
from matplotlib import pyplot as plt
import numpy.linalg as LA

iris = datasets.load_iris()
iris_df = pd.DataFrame(iris["data"], columns = iris["feature_names"])
iris_df["target"] = iris["target"]
group_df = [] #種類ごとにグループに分ける
N = [] #各グループのデータの数
mu = [] #各グループの標本平均
cov = [] #各グループの標本共分散行列
G = 3 #種類の数
for i in range(G):
  group_df.append(iris_df[iris_df["target"] == i])
  N.append(group_df[i].shape[0])
  mu.append(np.mean(group_df[i].drop(["target"], axis = 1)))
  cov.append(np.cov(group_df[i].drop(["target"], axis = 1).T, bias = True))

all_mu = np.mean(iris_df.drop(["target"], axis = 1)) #全標本の標本平均
matrix_B = 0 #群間変動行列
matrix_W = 0 #群内変動行列
for i in range(G):
  matrix_B = matrix_B + N[i] * np.dot(np.matrix(mu[i] - all_mu).T, np.matrix(mu[i] - all_mu))
  matrix_W = matrix_W + (N[i] - 1) * cov[i]

#固有値問題として帰着
S = np.dot(np.linalg.inv(matrix_W), matrix_B)
w, v = LA.eig(S)
sorted_index = np.argsort(w)
v_1, v_2 = np.array((v[sorted_index])[-1]), np.array((v[sorted_index])[-2])

df_x = np.dot(iris_df.drop(["target"], axis = 1), v_1.T)
df_y = np.dot(iris_df.drop(["target"], axis = 1), v_2.T)
color = []
color.append(np.repeat(["red"], N[0]))
color.append(np.repeat(["blue"], N[1]))
color.append(np.repeat(["green"], N[2]))
color = np.concatenate((color[0], color[1], color[2]))
plt.scatter(df_x, df_y, c = np.array(color))

f:id:mt19:20181021142101p:plain
それなりに可視化できている気がします。別に2次元じゃなくても良いよねと思ったので、次は3次元で可視化してみます。

from mpl_toolkits.mplot3d import Axes3D
v_3 = np.array((v[sorted_index])[-3])
df_z = np.dot(iris_df.drop(["target"], axis = 1), v_3.T)
fig = plt.figure()
ax = Axes3D(fig)
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.scatter3D(df_x, df_y, df_z, c = color)
ax.view_init(15, 55)

f:id:mt19:20181021143208p:plain
3次元にするとより分離ができそうです。

フィッシャーの線形判別とマハラノビス距離による分類

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

今回はフィッシャーの線形判別を行ないます。簡単に言うと、2グループに分かれているデータを最も良く分離できそうな軸に射影して判別しましょうというものです。最も良い軸とは、2グループの標本平均差が大きく、かつ2グループのグループ内分散が小さいものとなります。もっと視覚的に言うと、射影したらそれぞれのグループがこじんまりしててしかもお互いが離れるような軸を探せということになります。理論的な詳細はいくらで出てくると思うので、実装します。

まずは定番ロジスティック回帰でしてみます。データは昔作ったのと同じです。

from sklearn import linear_model
import numpy as np
import pandas as pd
from matplotlib import pyplot as plt

#create data
x_1 = 5
y_1 = 7
x_2 = 7
y_2 = 5
x1 = x_1 + np.random.randn(50)
y1 = y_1 + np.random.randn(50)
x2 = x_2 + np.random.randn(50)
y2 = y_2 + np.random.randn(50)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2]))
y = pd.Series(np.hstack([y1, y2]))
c = pd.Series(np.repeat([1, 0], 50))
data = pd.DataFrame({"x":x, "y":y, "c":c})

#logistic regression
reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga").fit(data[["x", "y"]], data["c"])
print(reg.intercept_, reg.coef_[0]) #->[-0.05453222] [-1.39388269  1.36970588]

#plot decision boundary
seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1)
plt.scatter(data["x"][0:49], data["y"][0:49], c = "red")
plt.scatter(data["x"][50:], data["y"][50:], c = "blue")
plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_ + reg.coef_[0][0] * seq), c = "black")
plt.grid(True)
plt.xlabel("x")
plt.ylabel("y")
plt.title("logistic regression")

#accuracy
print(reg.score(data[["x", "y"]], data["c"])) #->0.9でした

f:id:mt19:20181019202719p:plain
100個データがあるので10個判定を間違えているようです。ただ直線を引いただけなのでこれくらいでしょう。次にフィッシャーの線形判別を行ないます。

#グループ0かどうかを判定する
def fisher_function(x, x_1to2, w):
  output = np.dot(w.T, x) - (1 / 2) * np.dot(w.T, x_1to2)
  if output >= 0:
    return 1
  else:
    return 0


#データの整備
n1 = 50
n2 = 50
x_mean_11 = np.mean(x[0:n1])
x_mean_12 = np.mean(y[0:n1])
x_mean_21 = np.mean(x[n1:n1 + n2])
x_mean_22 = np.mean(y[n1:n1 + n2])
x_1_2 = np.array([x_mean_11, x_mean_12]) - np.array([x_mean_21, x_mean_22])
s_1to11 = np.sum(np.power(x[0:n1] - x_mean_11, 2))
s_1to12 = np.dot(x[0:n1] - x_mean_11, y[0:n1] - x_mean_12)
s_1to22 = np.sum(np.power(y[0:n1] - x_mean_12, 2))
s_2to11 = np.sum(np.power(x[n1:n1 + n2] - x_mean_21, 2))
s_2to12 = np.dot(x[n1:n1 + n2] - x_mean_21, y[n1:n1 + n2] - x_mean_22)
s_2to22 = np.sum(np.power(y[n1:n1 + n2] - x_mean_22, 2))
S_1 = np.array([[s_1to11, s_1to12], [s_1to12, s_1to22]])
S_2 = np.array([[s_2to11, s_2to12], [s_2to12, s_2to22]])
S = ((n1 - 1) * S_1 + (n2 - 1) * S_2) / (n1 + n2 - 2)

#射影する軸の方向ベクトル
w_hat = np.dot(np.linalg.inv(S), x_1_2)

#accuracyを求めてみる
x_1to2 = np.array([x_mean_11, x_mean_12]) + np.array([x_mean_21, x_mean_22])
input_data = np.array([x, y])
output_data = []
for i,j in zip(input_data[0], input_data[1]):
  output_data.append(fisher_function([i,j], x_1to2, w_hat))
print(np.sum(c == output_data) / (n1 + n2)) #->0.9でした

#最後に可視化
y_line = w_hat[0] * x + w_hat[1] * y
x_line = np.concatenate((y_line[0:n1] * w_hat[0] / w_hat[1],y_line[n1:n1 + n2] * w_hat[0] / w_hat[1]))
plt.scatter(x_line[0:n1], y_line[0:n1], c = "red")
plt.scatter(x_line[n1:n1 + n2], y_line[n1:n1 + n2], c = "black")
seq = np.arange(np.min(x_line) - 0.5, np.max(x_line) + 0.5)
plt.plot(seq, seq * w_hat[1] / w_hat[0], c = "blue")
plt.ylim(np.min(y_line) - 0.1, np.max(y_line) + 0.1)
plt.xlim(np.min(x_line) - 0.1, np.max(x_line) + 0.1)
plt.xlabel("x")
plt.ylabel("y")
plt.grid(True)

f:id:mt19:20181019202947p:plain
たまたまだと思いますがロジスティック回帰と同じaccuracyになりました。このように軸に射影して真ん中の点で分類といういたってシンプルなのが特徴です。行列を駆使すればもっとシンプルに書けると思います。すみませんgdgdで。

次に多クラスの分類をしてみます。これについてもフィッシャーの線形判別でできるそうなのですが、今回はマハラノビス距離による分類を行ないます。ただ類似度としてマハラノビス距離を利用しただけです。さすがに上のコードが自分で見てもひどいなと思ったので今回は短く抑えました。

#create data
n = 50
x_1 = 0
y_1 = 5
x_2 = 6
y_2 = 1
x_3 = 9
y_3 = 10
x1 = x_1 + np.random.randn(n)
y1 = y_1 + np.random.randn(n)
x2 = x_2 + np.random.randn(n)
y2 = y_2 + np.random.randn(n)
x3 = x_3 + np.random.randn(n)
y3 = y_3 + np.random.randn(n)
#上を元にデータをpandas.dataframeで作成
x = pd.Series(np.hstack([x1, x2, x3]))
y = pd.Series(np.hstack([y1, y2, y3]))
c = pd.Series(np.repeat([2, 1, 0], n))
data = pd.DataFrame({"x":x, "y":y, "c":c})
plt.scatter(data["x"][0:n], data["y"][0:n], c = "red")
plt.scatter(data["x"][n:2 * n], data["y"][n:2 * n], c = "blue")
plt.scatter(data["x"][2 * n:], data["y"][2 *n:], c = "green")

#logistic regression
reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga").fit(data[["x", "y"]], data["c"])
print(reg.intercept_)
print(reg.coef_)

#plot decision boundary
seq = np.arange(data["x"].min() - 1, data["x"].max() + 1, 0.1)
plt.scatter(data["x"][0:n], data["y"][0:n], c = "red")
plt.scatter(data["x"][n:2 * n], data["y"][n:2 * n], c = "blue")
plt.scatter(data["x"][2 * n:], data["y"][2 *n:], c = "green")
plt.plot(seq, (-1 / reg.coef_[0][1]) * (reg.intercept_[0] + reg.coef_[0][0] * seq), c = "black")
plt.plot(seq, (-1 / reg.coef_[1][1]) * (reg.intercept_[1] + reg.coef_[1][0] * seq), c = "purple")
plt.plot(seq, (-1 / reg.coef_[2][1]) * (reg.intercept_[2] + reg.coef_[2][0] * seq), c = "pink")
plt.xlabel("x")
plt.ylabel("y")
plt.title("logistic regression")
plt.grid(True)

#accuracy
print(reg.score(data[["x", "y"]], data["c"])) #->1.0

f:id:mt19:20181019220343p:plain
ロジスティック回帰でやるとこんな感じ。これをマハラノビス距離で行ないます。

#マハラノビス距離による多クラス分類
group = []
mean = []
cov = []
num_of_group = 3
for i in range(num_of_group):
  group.append(np.array(data[["x", "y"]][i * n:(i + 1) * n]).T)
  mean.append(np.mean(group[i], axis = 1))
  cov.append(np.cov(group[i], bias = True))
S = (1 / (n * num_of_group - num_of_group)) * (n - 1) * (cov[0] + cov[1] + cov[2])

#最短のマハラノビス距離となるグループを求める
D_g = []
S_inv = np.linalg.inv(S)
for g in range(num_of_group):
  for i in range(n):
    D_list = []
    for k in range(num_of_group):
      x = (group[2 - g].T)[i] - mean[k]
      D_list.append(np.dot(np.dot(x.T, S_inv), x))
    D_g.append(D_list.index(np.min(D_list)))
print(np.sum(np.array(D_g) == data["c"]) / len(D_g))#->1.0

フィッシャーの線形判別と比べてスッキリさせました。グループがはっきりしているのでもちろん100%分類できています。2次元から多次元にするのも容易です。

XORのようなデータを非線形ロジスティック回帰で分類する

参考書は相変わらずこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

今回の手順は前回とそっくりなので前回のを見てもらうとだいたい分かります。
基底関数展開を使った非線形回帰 - lazy diary
前回は基底関数展開から非線形の単回帰を行ないました。今回は、基底関数展開を利用した非線形のロジスティック回帰を行ないます。基底関数を求めるところまでは同じで、最後がロジスティック回帰になっています。
非線形の恩恵を確かめるために、通常の線形ロジスティック回帰では分類できないXORのようなデータを今回は使います。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.cluster import KMeans
from sklearn import linear_model
from scipy.stats import multivariate_normal

cov = 0.5 * np.identity(2)
data1 = pd.DataFrame(multivariate_normal.rvs(mean = (-2.0, -2.0), cov = cov, size = 25))
data2 = pd.DataFrame(multivariate_normal.rvs(mean = (2.0, -2.0), cov = cov, size = 25))
data3 = pd.DataFrame(multivariate_normal.rvs(mean = (2.0, 2.0), cov = cov, size = 25))
data4 = pd.DataFrame(multivariate_normal.rvs(mean = (-2.0, 2.0), cov = cov, size = 25))

data1["c"] = 0
data2["c"] = 1
data3["c"] = 0
data4["c"] = 1

data = pd.concat([data1, data2, data3, data4], join = "inner")
data = data.reset_index(drop = True)
data = data.rename(columns = {0:"x", 1:"y", "c":"c"})

x = data["x"]
y = data["y"]
clusters = data["c"]

plt.scatter(x, y, c = clusters)

f:id:mt19:20181109161138p:plain
色付けをまともにしてないので見づらいですが、たしかにXORです。ここから前回と同じようにk-meansからの基底関数の構成およびデータの変換を行ないます。

#k-meansなどでクラスタリングする
dist = []
for i in range(1,30):
  km = KMeans(n_clusters = i)
  km.fit(data[["x", "y"]])
  dist.append(km.inertia_)
  
plt.plot(range(1,30), dist, marker = ".")

f:id:mt19:20181109161217p:plain
エルボー則に従って4くらいかなと考えられます。

from sklearn import linear_model
#クラスタ数を4とする

#ガウス型基底関数
def gauss_func(x, mu, sigma):
  return np.exp((np.sum(np.power(x - mu, 2))) / (-2 * sigma))

m = 4
km = KMeans(n_clusters = m)
km.fit(data[["x", "y"]])
cluster_center = km.cluster_centers_
dist = np.power(km.transform(data[["x", "y"]]), 2)
cluster_num = km.predict(data[["x", "y"]])
least_dist = np.zeros(shape = (m, 1))
for i in range(len(data)):
  for j in range(m):
    if cluster_num[i] == j:
      least_dist[j] += dist[i][j]
for i in range(m):
  least_dist[i] = least_dist[i] / np.sum(cluster_num == i)

transformed_X = pd.DataFrame(np.zeros(shape = (len(data), m)))
for i in range(len(data)):
  for j in range(m):
    transformed_X.iloc[i, j] = gauss_func(data[["x", "y"]].iloc[i,:], cluster_center[j], least_dist[j])

さて、カーネルパラメータの推定とカーネルを使ってデータの変換まで終わったのでロジスティック回帰を実行します。

#ロジスティック回帰の実行
log_reg = linear_model.LogisticRegression(penalty = "l2", solver = "saga", C = 1000.0, max_iter = 1000).fit(transformed_X, clusters)
print(log_reg.score(transformed_X, clusters)) #->1.0

とりあえず分類できました。決定境界をプロットしてみます。

seq = np.linspace(data["x"].min(), data["x"].max(), 50)
X1, X2 = np.meshgrid(seq, seq)
Z = []
x1 = X1.reshape(-1, 1)
x2 = X2.reshape(-1, 1)
testdata = np.hstack((x1, x2))
transformed_testdata = pd.DataFrame(np.zeros(shape = (len(x1), m)))
for i in range(len(x1)):
  for j in range(m):
    transformed_testdata.iloc[i, j] = gauss_func(testdata[i,:], cluster_center[j], least_dist[j])
Z = log_reg.predict_proba(transformed_testdata)

plt.title("Non Lnear Soft Margin SVM")
plt.xlabel("x")
plt.ylabel("y")
Z = np.array(Z[:,0]).reshape(X1.shape)
plt.contour(X1, X2, Z, colors = "black", levels = [0.5])
plt.scatter(data["x"], data["y"], c = clusters)

f:id:mt19:20181109161626p:plain
ガウスカーネルっぽい丸みを帯びた決定境界です。

基底関数展開を使った非線形回帰

参考書はこちら

多変量解析入門――線形から非線形へ

多変量解析入門――線形から非線形へ

この第3章非線形回帰モデルの内容です。

今回は基底関数展開というものを行ないます。基底というのは線形代数で学ぶ基底と同じ意味合いです。ここでは基底関数ですので、その基底の要素が関数になっているわけです。例えば非線形回帰モデルで最も有名な例として、d次元多項式回帰を考えると、その基底関数は\{1,x,x^2,x^3,\dots,x^d\}となります。この関数を適用したデータに対して重回帰モデルをfitさせれば非線形モデリングを行なうことが出来ます。
今回は、多項式回帰ではなく基底関数としてガウス型基底関数というものを用います。これはどんな関数かというとガウス分布の密度関数のカーネル部分です。方法としてはざっくり次のようになります。

1. 基底関数の個数を決めるためにデータに対してk-meansなどでクラスタリングを行なう。
2 そのクラス数を基底関数の個数とし、また各クラスに属するデータの平均と距離による分散を求め、それをそのクラスのガウス型基底関数のパラメータとする。
3. 各データにガウス型基底関数をそれぞれapplyし、そのデータに対して重回帰を行なう。

とりあえずやってみます。データはいかにも線形で表せないようなデータを自分で作ります。

import numpy as np
import pandas as pd
from matplotlib import pyplot as plt
from sklearn.cluster import KMeans
from sklearn import linear_model
from scipy.stats import norm
#データのもととなる非線形関数
def f(x):
  return ((x**3 + 2 * (x**2) - 24 * x) / 10) + np.sin(2 * x)

#可視化
seq = np.arange(-5, 5, 0.2)
plt.scatter(seq, f(seq), c = "red")

X = seq + norm.rvs(loc = 0, scale = 1, size = len(seq))
y = f(X) + norm.rvs(loc = 0, scale = 1, size = len(seq))
plt.scatter(X, y, c = "blue")

f:id:mt19:20181011101211p:plain
赤色が今回使う非線形関数です。そこにノイズを加えたのが青色のデータです。問題としては青色しか分かっていない状態が想定されています。
さて、ここからk-meansを当てはめます。

X = X.reshape((-1, 1))
#k-meansなどでクラスタリングする
dist = []
for i in range(1,8):
  km = KMeans(n_clusters = i)
  km.fit(X)
  dist.append(km.inertia_)
  
plt.plot(range(1,8), dist, marker = ".")

f:id:mt19:20181011101446p:plain
クラス数はエルボー則に基づき、2としても良いのですがどこか心もとないので3にします。

#クラスタ数は3とする
m = 3
km = KMeans(n_clusters = m)
km.fit(X)
cluster_center = km.cluster_centers_ #平均
dist = np.power(km.transform(X), 2)
cluster_num = km.predict(X)
h = [0, 0, 0] #分散
for i in range(len(X)):
  if cluster_num[i] == 0:
    h[0] += dist[i][0]
  elif cluster_num[i] == 1:
    h[1] += dist[i][1]
  else:
    h[2] += dist[i][2]
h[0] = h[0] / np.sum(cluster_num == 0)
h[1] = h[1] / np.sum(cluster_num == 1)
h[2] = h[2] / np.sum(cluster_num == 2)

これで、3つのガウス型基底関数のパラメータが求められました。さて、元のデータにこの関数をそれぞれapplyします。

#ガウス型基底関数
def basis_function(x, mu, h):
  return np.exp(np.power(x - mu, 2) / (-2 * h))


#値を変換する
transformed_X = pd.DataFrame(np.zeros(shape = (len(seq), m)))
for i in range(m):
  transformed_X.iloc[:,i] = basis_function(X, cluster_center[i], h[i])

これで準備完了です。最後に重回帰させ、結果をプロットします。

#重回帰させる
lm = linear_model.LinearRegression(fit_intercept = True, normalize = False)
lm.fit(transformed_X, y)
plt.scatter(X, lm.predict(transformed_X), c = "red")
plt.scatter(X, y, c = "blue")

f:id:mt19:20181011101857p:plain
赤色が今回得られた推定値です。青色がfitさせたもとのデータです。なんとか3つの基底関数で表現するとこのくらいになるのかなといった感じです。
余力があればそのうち説明変数を増やしたバージョンもやるかもしれません。

追記1:
私がひねくれていたのが良くなかったか、素直にエルボー則に従ってクラスタ数2で実行したところ次のようになりました。
f:id:mt19:20181011103054p:plain

追記2:
たまたまscikit-learnのチートシートを眺めていたらこれを簡単に実装できそうなメソッドがありそうなのでやってみます。それはsupport vector regressionというものです。数理的な内容はここが簡潔にまとまっています。
https://datachemeng.com/supportvectorregression/
最適化の部分でL2正則化項が入っていたり、誤差関数が少し違ったりと色々と細かい部分で違いはあるのですが、これらハイパーパラメータをいじれば基底関数展開に近づけることができるのでさほど問題はありません。やってみましょう。(データは上と同じものを使用)

from sklearn import svm
svr = svm.SVR(kernel = "rbf",C = 1000, epsilon = 0.0001).fit(X, y)
plt.scatter(X, svr.predict(X), c = "red")
plt.scatter(X, y, c = "blue")

f:id:mt19:20181013101442p:plain
赤色が今回の予測値です。青色は元のデータです。上の基底関数展開に近づけるため正則化項はほぼ無効化させています。これだけで実装できるのだから便利ですね。