はじめに

以前、「機械学習の分類」で取り上げたアルゴリズムについて、その理論とpythonでの実装、scikit-learnを使った分析についてステップバイステップで学習していく。個人の学習用として書いてるので間違いなんかは大目に見て欲しいと思います。

今回はロジスティック回帰について。ロジスティック回帰も、回帰と書いてはいるものの、パーセプトロンのように二値分類を扱うアルゴリズムです。

今回参考にしたのは以下のサイト。ありがとうございます。

理論

ロジスティック回帰の理論について、まずは活性化関数であるシグモイド関数を導出してみます。

シグモイド関数

ロジスティック回帰は二値分類なので、クラス$C_1$と$C_2$について考える。$C_1$の確率$P(C_1)$と$C_2$の確率$P(C_2)$の合計は1です。

データ列$\boldsymbol{x}$を与えた時に$C_1$になる確率は、ベイズの定理より

\begin{align}
P(C_1|\boldsymbol{x})&=\frac{P(\boldsymbol{x}|C_1)P(C_1)}{P(\boldsymbol{x})} \\
&= \frac{P(\boldsymbol{x}|C_1)P(C_1)}{P(\boldsymbol{x}|C_1)P(C_1)+P(\boldsymbol{x}|C_2)P(C_2)} \\
&= \frac{1}{1+\frac{P(\boldsymbol{x}|C_2)P(C_2)}{P(\boldsymbol{x}|C_1)P(C_1)}} \\
&= \frac{1}{1+\exp(-\ln\frac{P(\boldsymbol{x}|C_1)P(C_1)}{P(\boldsymbol{x}|C_2)P(C_2)})} \\
&= \frac{1}{1+\exp(-a)} = \sigma(a)
\end{align}

この$\sigma(a)$のことをシグモイド関数と呼びます。シグモイド関数は以下のように、0から1の値をとるので、確率を表すに都合の良い関数です。

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

与えられたデータ列$\boldsymbol{x}=(x_0,x_1,\cdots,x_n)$と教師のクラス分類$\boldsymbol{t}=(t_0,t_1,\cdots,t_n)$を用い、

L(\boldsymbol{x})=\frac{1}{1+\exp(-\boldsymbol{w}^T\boldsymbol{x})}

のパラメータ$\boldsymbol{w}=(w_0,w_1,\cdots,w_n)$を最適化していきます。

クロスエントロピー誤差

ある$x_i$が与えられたときにクラスが$C_1$になる確率$P(C_1|x_i)$を$p_i$とすると、クラスが$C_2$になる確率$P(C_2|x_i)$は$(1-p_i)$となる。つまり、クラスが$t_i$になる確率$P(t_i|x_i)$は、$$P(t_i|x_i)=p_i^{t_i}(1-p_i)^{1-t_i}$$となる。

これを全データに適用すると、

\begin{align}
P(\boldsymbol{t}|\boldsymbol{x})&=P(t_0|x_0)P(t_1|X_1)\cdots P(t_{n-1}|x_{n-1}) \\
&=\prod_{i=0}^{n-1}P(t_i|x_i) \\
&=\prod_{i=1}^{n-1}p_i^{t_i}(1-p_i)^{1-t_i}
\end{align}

となります。この両辺の対数をとると、

\log P(\boldsymbol{t}|\boldsymbol{x}) = \sum_{i=0}^{n-1}\{t_i\log p_i+(1-t_i)\log (1-p_i)\}

これは対数尤度と呼ばれ、対数尤度を最大にするために、符号を反転して、

E(\boldsymbol{x}) = -\frac{1}{n}\log P(\boldsymbol{t}|\boldsymbol{x}) = \frac{1}{n}\sum_{i=0}^{n-1}\{-t_i\log p_i-(1-t_i)\log (1-p_i)\}

この$E$をクロスエントロピー誤差関数と言います。
あとで使うので、$E$の微分は、

\frac{\partial{E}}{\partial{w_i}}=\frac{1}{n}\sum_{i=0}^{n-1}(p_i-t_i)x_i

となります(説明省略)

共役勾配法

さて、クロスエントロピー誤差関数を最小化するためには、以前も出てきた勾配法を使います。ここでも、最急降下法や確率的勾配降下法が使えるんですが、**共役勾配法(Conjugate Gradient Method)**を使います。詳しい話はWikipedia:共役勾配法に譲りますが、最急勾配法と比べると高速かつ学習率を設定しなくても収束するというアルゴリズムです。
これもpythonで実装したいところですが、面倒なので(!)pythonのscipy.optimize.fmin_cgというライブラリを使います。

pythonによる実装

ここまでの理論を使ってLogisticRegressionクラスを実装していきます。fmin_cgは勾配関数を与えると良好な結果が出ると言うことで、先ほどの勾配関数を使います。

from scipy import optimize

class LogisticRegression:
  def __init__(self):
    self.w = np.array([])

  def sigmoid(self, a):
    return 1.0 / (1 + np.exp(-a))
  
  def cross_entropy_loss(self, w, *args):
    def safe_log(x, minval=0.0000000001):
      return np.log(x.clip(min=minval))
    t, x = args
    loss = 0
    for i in range(len(t)):
      ti = (t[i]+1)/2
      h = self.sigmoid(w.T @ x[i])
      loss += -ti*safe_log(h) - (1-ti)*safe_log(1-h)

    return loss/len(t)

  def grad_cross_entropy_loss(self, w, *args):
    t, x = args
    grad = np.zeros_like(w)
    for i in range(len(t)):
      ti = (t[i]+1)/2
      h = self.sigmoid(w.T @ x[i])
      grad += (h - ti) * x[i]

    return grad/len(t)

  def fit(self, x, y):
    w0 = np.ones(len(x[0])+1)
    x = np.hstack([np.ones((len(x),1)), x])

    self.w = optimize.fmin_cg(self.cross_entropy_loss, w0, fprime=self.grad_cross_entropy_loss, args=(y, x))

  @property
  def w_(self):
    return self.w

このクラスを使ってirisのデータを分類してみます。境界もあわせて描画します。境界は$\boldsymbol{w}^T\boldsymbol{x}=0$の線です。2クラスを1と-1にしたのでそれに合うようなコードにしています。

df = df_iris[df_iris['target']!='setosa']
df = df.drop(df.columns[[1,2]], axis=1)
df['target'] = df['target'].map({'versicolor':1, 'virginica':-1})

# グラフの描画
fig, ax = plt.subplots()

df_versicolor = df_iris[df_iris['target']=='versicolor']

x1 = df_iris[df_iris['target']=='versicolor'].iloc[:,3].values
y1 = df_iris[df_iris['target']=='versicolor'].iloc[:,0].values

x2 = df_iris[df_iris['target']=='virginica'].iloc[:,3].values
y2 = df_iris[df_iris['target']=='virginica'].iloc[:,0].values

xs = StandardScaler()
ys = StandardScaler()

xs.fit(np.append(x1,x2).reshape(-1, 1))
ys.fit(np.append(y1,y2).reshape(-1, 1))

x1s = xs.transform(x1.reshape(-1, 1))
x2s = xs.transform(x2.reshape(-1, 1))
y1s = ys.transform(y1.reshape(-1, 1))
y2s = ys.transform(y2.reshape(-1, 1))

x = np.concatenate([np.concatenate([x1s, y1s], axis=1), np.concatenate([x2s, y2s], axis=1)])

y = df['target'].values

model = LogisticRegression()
model.fit(x, y)

ax.scatter(x1s, y1s, color='red', marker='o', label='versicolor')
ax.scatter(x2s, y2s, color='blue', marker='s', label='virginica')

ax.set_xlabel("petal width (cm)")
ax.set_ylabel("sepal length (cm)")

# 分類境界を描画する
w = model.w_
x_fig = np.linspace(-2.,2.,100)
y_fig = [-w[1]/w[2]*xi-w[0]/w[2] for xi in x_fig]
ax.plot(x_fig, y_fig)
ax.set_ylim(-2.5,2.5)

ax.legend()
print(w)
plt.show()

Optimization terminated successfully.
         Current function value: 0.166434
         Iterations: 12
         Function evaluations: 41
         Gradient evaluations: 41
[-0.57247091 -5.42865492 -0.20202263]

結構綺麗に分類できているようです。

scikit-learnの実装

scikit-learnもLogisticRegressionクラスがあるので、上のコードとほとんど同じです。

from sklearn.linear_model import LogisticRegression

df = df_iris[df_iris['target']!='setosa']
df = df.drop(df.columns[[1,2]], axis=1)
df['target'] = df['target'].map({'versicolor':1, 'virginica':-1})

# グラフの描画
fig, ax = plt.subplots()

df_versicolor = df_iris[df_iris['target']=='versicolor']

x1 = df_iris[df_iris['target']=='versicolor'].iloc[:,3].values
y1 = df_iris[df_iris['target']=='versicolor'].iloc[:,0].values

x2 = df_iris[df_iris['target']=='virginica'].iloc[:,3].values
y2 = df_iris[df_iris['target']=='virginica'].iloc[:,0].values

xs = StandardScaler()
ys = StandardScaler()

xs.fit(np.append(x1,x2).reshape(-1, 1))
ys.fit(np.append(y1,y2).reshape(-1, 1))

x1s = xs.transform(x1.reshape(-1, 1))
x2s = xs.transform(x2.reshape(-1, 1))
y1s = ys.transform(y1.reshape(-1, 1))
y2s = ys.transform(y2.reshape(-1, 1))

x = np.concatenate([np.concatenate([x1s, y1s], axis=1), np.concatenate([x2s, y2s], axis=1)])

y = df['target'].values

model = LogisticRegression(C=100)
model.fit(x, y)

ax.scatter(x1s, y1s, color='red', marker='o', label='versicolor')
ax.scatter(x2s, y2s, color='blue', marker='s', label='virginica')

ax.set_xlabel("petal width (cm)")
ax.set_ylabel("sepal length (cm)")

# 分類境界を描画する
w = model.coef_[0]

x_fig = np.linspace(-2.,2.,100)
y_fig = [-w[0]/w[1]*xi-model.intercept_/w[1] for xi in x_fig]
ax.plot(x_fig, y_fig)
ax.set_ylim(-2.5,2.5)

ax.legend()
plt.show()

こちらもうまい具合に分類できています。

まとめ

機械学習の世界でも割と重要（と思われる）なロジスティック回帰についてまとめました。だんだんこの辺りから理論が難しくなってきましたね。