入力を非線形変換したロジスティック回帰 - 機械学習に詳しくなりたいブログ

入力を非線形変換したパーセプトロンでも基底変換を実験しましたが、ロジスティック回帰でも試してみます。パーセプトロンで試したときは、変換 $\phi(\mathbf{x})$ は適当に決めてしまいましたが、今回はガウス基底を使ってみます。名前がついてるってだけで、適当に選んだことには違いないですが。

まず基本に立ち返って、基底関数を使って入力を非線形変換すると何が嬉しいのかを確認してみます。書籍「パターン認識と機械学習（上）」にわかりやすい図が載っていましたので、図を再現するようなコードを書いてみました。

f:id:opabinia2:20180819213954p:plain 左のグラフが、元のデータ点 $(x_1,x_2)$ の分布です。これをガウス基底

$\begin{eqnarray*} \displaystyle \phi_1 &=& \exp\left(-\frac{(x_1 - \mu_{11})^{2} + (x_2 - \mu_{12})^{2}}{2\sigma^{2}}\right) \tag{1} \\ \displaystyle \phi_2 &=& \exp\left(-\frac{(x_1 - \mu_{21})^{2} + (x_2 - \mu_{22})^{2}}{2\sigma^{2}}\right) \tag{2} \\ \end{eqnarray*}$

を使って変換したものが右のグラフです。基底の中心は $(\mu_{11},\mu_{12})=(-1,-1)$ 、 $(\mu_{21},\mu_{22})=(0,0)$ としています。左のグラフの緑の線が基底の等高線です。基底の中心からの距離（のようなもの）が新たな軸 $\phi$ となりますので、元のデータでは離れた位置に分布している青色の点も、 $\phi_2$ の軸で見れば、中心からの距離が同じなので、基底変換後は縦軸ではかたまった分布となっています。一方 $\phi_1$ の軸で見れば、距離が離れていますから基底変換後の横軸では離れた位置に分布しています。このような変換を行うことで、 $\mathbf{w}^{T}\mathbf{x}$ では線形分離できなかった分布が、 $\mathbf{w}^{T}\mathbf{\phi(\mathbf{x})}$ では線形分離することができるようになりました。複雑な分布のデータも、適切な変換を行うことで線形識別のモデルが適用できるということですね。

ということで、以下のような線形で分離できないようなデータもロジスティック回帰で識別できることが確認できました。 f:id:opabinia2:20180819215909p:plain 基底は16個使っていて、中心や $\sigma$ は適当に決めたんですが、この取り方によっては計算がうまくいかないことがありました。調べてみると、シグモイド関数のオーバーフローが原因のようでした。プログラムの動作を追って原因を見つけたのではなく、webで調べて解決したんですが。この辺りはいろいろ対策や工夫があるようです。まあまずはこういう現象を実際に体験できたことが重要ですね。計算がうまくいかないとき、シグモイド関数への入力を調べみるっていう発想ができるので。というか数値計算では常識レベルの話なのかな。

関連：ロジスティック回帰を最急降下法とニュートン法で解く

今回のコードです。

# 入力を非線形変換したロジスティック回帰
import matplotlib.pyplot as plt
import numpy as np
from matplotlib import cm


# ガウス基底
def gaussian_basis(x1, x2, mu1, mu2):
    s = 1.0
    return np.exp(-((x1 - mu1)**2 + (x2 - mu2)**2)/(2.0*(s**2)))


# 基底変換
def phi(x1, x2):
    return np.array([gaussian_basis(x1, x2, mu0, mu1) for (mu0, mu1) in center])


# 判別結果
def f(x1, x2):
    if sigmoid(np.dot(w.T, phi(x1, x2))) < 0.5:
        return 0
    else:
        return 1


# 交差エントロピー誤差
def Ew(y, t):
    ew = 0
    for i in range(N):
        print(y[i])
        ew -= np.log(t[i]*y[i] + (1-t[i])*(1-y[i]))
    return ew


# 教師データの属するクラスを返す
def teaching(x, y):
    if ((x-0.5)**2 + (y-0.5)**2 < 0.25**2 and x < 0.5) or ((x-0.5)**2 + (y-0.5)**2 > 0.25**2 and x > 0.5) :
        return 0
    else:
        return 1


# シグモイド関数
def sigmoid(a):
    return 1.0/(1 + np.exp(-a))


# 重み付け行列Rを計算
def R_(y):
    R_ij = np.empty(N)
    for i in range(N):
        R_ij[i] = y[i] * (1 - y[i])
    return np.diag(R_ij)


# ランダムシードを固定
np.random.seed(0)

# 訓練データ数
N = 200

# 入力次元数M
M0 = 4
M = M0**2

# クラス数
K = 2

# 教師データ作成
x = np.random.uniform(0, 1, [N, 2])
t = np.vectorize(teaching)(x[:, 0], x[:, 1])
# 1 of K符号
T = np.eye(K)[t]

plt.xlim(0, 1)
plt.ylim(0, 1)
plt.scatter(x[t == 0, 0], x[t == 0, 1], color="blue", alpha=0.5)
plt.scatter(x[t == 1, 0], x[t == 1, 1], color="red", alpha=0.5)
plt.show()

# ガウス基底の中心点
center = np.array([[x0, x1] for x0 in range(M0) for x1 in range(M0)])

# 行列X(ガウス基底を通す）
X = np.empty([N, M])
for i in range(N):
    X[i, :] = phi(x[i, 0], x[i, 1])


# パラメータw初期化
w = np.random.uniform(-1, 1, M)

# 出力y
y = np.vectorize(sigmoid)(np.dot(w, X.T))

# 重み付け行列R
R = R_(y)

# ヘッセ行列
H = np.dot(np.dot(X.T, R), X)

# ニュートン法
for i in range(7):
    w -= np.dot(np.dot(np.linalg.inv(H), X.T), (y -t))
    # yを再計算
    y = np.vectorize(sigmoid)(np.dot(w, X.T))
    # R,Hを再計算
    R = R_(y)
    H = np.dot(np.dot(X.T, R), X)


# グラフ表示用の判別結果
a, b = np.meshgrid(np.linspace(0, 1, 500), np.linspace(0, 1, 500))
vec_f = np.vectorize(f)
c = vec_f(a, b)

plt.xlim(0, 1)
plt.ylim(0, 1)
plt.scatter(x[t == 0, 0], x[t == 0, 1], color="blue", alpha=0.5)
plt.scatter(x[t == 1, 0], x[t == 1, 1], color="red", alpha=0.5)
plt.contourf(a, b, c, alpha=0.2, cmap=cm.coolwarm)
plt.show()