機械学習のアルゴリズム（ロジスティック回帰）

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

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

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

* [NumPyでロジスティック回帰を実装する](https://ohke.hateblo.jp/entry/2019/06/29/230000)
* [ロジスティック回帰](https://masamunetogetoge.com/classification-logisticregression)
* [ロジスティック回帰式の導出](https://qiita.com/hatnychan/items/b39f8b166b5c5ec86c0e)
* [共役勾配法によるロジスティック回帰のパラメータ推定](http://aidiary.hatenablog.com/entry/20140415/1397570262)
* [交差エントロピー誤差関数](https://python.atelierkobato.com/cce/)

# 理論
ロジスティック回帰の理論について、まずは活性化関数である**シグモイド関数**を導出してみます。
## シグモイド関数
ロジスティック回帰は二値分類なので、クラス$C_1$と$C_2$について考える。$C_1$の確率$P(C_1)$と$C_2$の確率$P(C_2)$の合計は1です。

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

```math
\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の値をとるので、確率を表すに都合の良い関数です。

![sigmoid](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/296113/1d00eeec-3842-fc89-e2da-badab729ff62.png)

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

```math
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}$$となる。

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

```math
\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}
```
となります。この両辺の対数をとると、

```math
\log P(\boldsymbol{t}|\boldsymbol{x}) = \sum_{i=0}^{n-1}\{t_i\log p_i+(1-t_i)\log (1-p_i)\}
```
これは**対数尤度**と呼ばれ、対数尤度を最大にするために、符号を反転して、

```math
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$の微分は、

```math
\frac{\partial{E}}{\partial{w_i}}=\frac{1}{n}\sum_{i=0}^{n-1}(p_i-t_i)x_i
```
となります(説明省略)

## 共役勾配法
さて、クロスエントロピー誤差関数を最小化するためには、以前も出てきた勾配法を使います。ここでも、最急降下法や確率的勾配降下法が使えるんですが、**共役勾配法(Conjugate Gradient Method)**を使います。詳しい話は[Wikipedia:共役勾配法](https://ja.wikipedia.org/wiki/%E5%85%B1%E5%BD%B9%E5%8B%BE%E9%85%8D%E6%B3%95)に譲りますが、最急勾配法と比べると高速かつ学習率を設定しなくても収束するというアルゴリズムです。
これもpythonで実装したいところですが、面倒なので(!)pythonの[scipy.optimize.fmin_cg](https://docs.scipy.org/doc/scipy-0.13.0/reference/generated/scipy.optimize.fmin_cg.html)というライブラリを使います。

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

```python
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にしたのでそれに合うようなコードにしています。

```python
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]
```

<img width="388" alt="logistic_regression_1.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/18497/2f8f2fd7-b1ad-5d3a-7da6-12e35a1d3c46.png">

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

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

```python
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()
```

<img width="390" alt="logistic_regression_2.png" src="https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/18497/a4478b00-8470-c9a1-c666-9d18353b6fdd.png">

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

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