ロジスティック回帰を実装してみよう

こんにちは、小澤です。

当エントリは「Machine Learning Advent Calendar 2017」の2日目のエントリです。

今回は、もっとも基本的な分類の1つであるロジスティック回帰を実装してみます。

分類の話

分類問題では、データを分類するための線(決定境界という)ものを決めます。 これはどういうことかというと、正解が「◯」か「×」かの二値分類を考えてみましょう。

以下のデータはx軸とy軸の2つの特徴量(テーブルデータ中の列)で表されるとします。

このデータに対して、例えば以下のように線を引きます。

この線より下側は1, 上側は0を出力するような関数を作って、0なら「◯」、1なら「×」というような出力を行うのが分類問題になります。 分類問題では、この線をどのように引くかが手法ごとの特徴となります。

手法によっては以下のデータのような直線ではあらわせないものを曲線などで分けるものあります。

また、実際のデータでは、以下のような完全に分離できない状況のものがほとんどです。

これを曲線で分離できる手法を使って無理やり全部分類できるようにしようとすると、学習時に利用したデータしか使えない「過学習」という状況に陥ります。

機械学習では、すでに手に入っているデータで学習を行って作成したモデルを、これから入力される未知のデータに適用することでその結果を予測するということが行われます。 そのため、過学習していないモデルを作成することも重要になります。

ロジスティック回帰について

ロジスティック回帰がどのようなものなのかは、以前解説しているので、そちらも合わせてご覧ください。

• ロジスティック回帰をしてみる #alteryx #11 | Alteryx Advent Calendar 2016 ｜ Developers.IO

実装してみる

さて、実装してみましょう。

利用するデータ

今回利用するデータは、前回の冒頭でも紹介したirisとなります。 また、ロジスティック回帰は二値分類となりますが、データセット中の品種は3種類ありますので、そのうち2種類のみを利用します(one-versus-rest法や多項ロジットモデルは利用しません)。

まずは、データを確認してみます。

from sklearn.datasets import load_iris
iris = load_iris()
iris_df = pd.DataFrame(iris.data, columns=iris.feature_names)
iris_df['label'] = iris.target
iris_df

正解のラベルには品種ごとに0, 1, 2のいずれかが入っています。

前処理

まず前回同様切片の項を追加してみましょう

# 切片を追加
iris_data = np.zeros([iris.data.shape[0], iris.data.shape[1] + 1])
for i in range(iris.data.shape[0]):
    for j in range(iris.data.shape[1]):
        iris_data[i, j] = iris.data[i, j]
    iris_data[i, iris.data.shape[1]] = 1
iris_data = iris_data[50:]

続いて、今回用に二値分類にするために、データを一部のみにします。 このデータセットは50個ずつ並んでいるので、先頭50行を読み飛ばすことで後半2つの二値分類にできます。

iris_data = iris_data[50:]
# 後半２つのラベルは1と2だが、これを0と1にしてる
iris_target = iris.target[50:] - 1

次に、学習データとテストデータにします。 先ほども書いた通り、このデータセットは正解ラベルごとに並んでいるので、分割に際して適当な位置で切ってしまうと、偏りが発生します。 今回は、一度ランダムにならべかえてから前半70個を学習用、後半30個をテスト用に分割しています。

# ランダムなインデックス作成
index = np.arange(len(iris_data))
np.random.shuffle(index)

# 学習データは前半70個
iris_train_data =  iris_data[index[:70]]
iris_train_target = iris_target[index[:70]]

# テストデータは後半30個
iris_test_data = iris_data[index[70:]]
iris_test_target = iris_target[index[70:]]

そのあとは、線形回帰の時と同様、データの正規化をしています。

# 各列の平均
iris_mu = [sum(iris_train_data[:, row]) / len(iris_train_data[:, row]) for row in range(iris_train_data.shape[1])]

# 各列の標準偏差
iris_sigma = [math.sqrt(sum((iris_train_data[:, row] - iris_mu[row]) ** 2  / len(iris_train_data[:, row]))) for row in range(iris_train_data.shape[1])]

# 正規化
for i in range(iris_train_data.shape[0]):
    for j in range(iris_train_data.shape[1] - 1):
        iris_train_data[i, j] = (iris_train_data[i, j] - iris_mu[j]) / iris_sigma[j]
        
for i in range(iris_test_data.shape[0]):
    for j in range(iris_test_data.shape[1] - 1):
        iris_test_data[i, j] = (iris_test_data[i, j] - iris_mu[j]) / iris_sigma[j]

学習

さて、いよいよ学習です。 とは言っても、ほとんどは線形回帰の時と同じです。

まず、weightをランダムな値で初期化します。

w = random(len(iris.feature_names) + 1)

続いて、現在のweightでの予測をします。

[latex] \hat{y} = \frac{1}{1 + \exp^{(-{\bf w}^T{\bf x})}} [/latex]

y_hat = 1 / (1 + np.exp(-np.dot(iris_train_data, w)))

この値を元に、損失関数の計算をするのですが、今回もループ回数を固定にするため、この計算はしていません。 また、正則化にはL2正則化を利用しています。

[latex] L = \frac{1}{m} (-\sum^{m}_{i} y_i \log(\sigma({\bf w}^T{\bf x}_i)) + (1 -y_i) \log(1-\sigma({\bf w}^T{\bf x}_i)) + \eta || {\bf w}||^2) \\ \frac{\partial E}{\partial {\bf w}_n} = -\frac{1}{m} (\sum^{m}_{i}\sigma(-{\bf w}^T{\bf x}){\bf x}_i + \eta {\bf w}_n) [/latex]

をイテレーションします。

for i in range(100):
    y_hat = 1 / (1 + np.exp(-np.dot(iris_train_data, w)))
    w = w +  alpha * ((np.dot((y_hat - iris_train_target), iris_train_data) - (eta * w)) / iris_train_data.shape[0])

テストデータで性能評価

最後にテストデータでの性能評価をしてみましょう。 正解率とPrecision, Recall, F値を求めています。

sum = 0
tp = tf = fn = fp = 0
for i, j in zip(iris_test_target, result):
    if i == j: sum += 1
    if i == 1 and j == 1:  tp += 1
    if i == 1 and j ==  0:  fn += 1
    if i == 0 and j == 1 : fp += 1
    if i == 0 and j == 0 : tn += 1
        
accuracy = sum / len(result)
precision  = tp / (tp + fp)
recall = tp / (tp + fn)
f_measure = (2 * precision * recall) / (precision + recall)
print('accuracy : ', accuracy)
print('precision : ', precision)
print('recall : ', recall)
print('f-measure : ', f_measure)

結果は以下のようになりました(初期値がランダムのため、異なる結果になる可能性もあります)。

accuracy :  0.8666666666666667
precision :  0.9230769230769231
recall :  0.8
f-measure :  0.8571428571428571

おわりに

今回は、ロジスティック回帰の実装をしてみました。

機械学習ではかなり簡単な部類に入るものですが、 一度やってみると、他の手法を実装するに際しても必要となる一連の流れが理解できるようなないようになっているかと思います。

明日は、主成分分析を実装してみます。 お楽しみに！