Python
numpy
MachineLearning
matplotlib

Courseraの機械学習のコースの宿題をPythonでやってみた②

前回の記事の続き.Courseraの機械学習のコースの宿題をTensorFlowとかは使わずにPythonとNumPyでトライします.

Prerequisite

  • Python 3.6.4
  • NumPy 1.13.3
  • Matplotlib 2.0.2

Programming Exercise 2: Logistic Regression

今回の問題は,あなたが大学の学部のadministratorであるとして,応募者が受験した2つの試験の点数を基に応募者の合格可能性をロジスティック回帰で求めようというものです

データとしては,過去の応募者の2つの試験の点数と合否が与えられています.ロジスティック回帰で応募者の合格可能性を計算します.

宿題を解く

今回も,各ポイントごとに解説を行なっていきます.

まず,必要になるライブラリ・モジュールをインポート:

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

データを読み込んで散布図を描く

まず,前回同様,numpy.loadtxtを用いてデータファイルを読み込み,2つの試験の点数データ$X$と合否データ$y$に分けます:

data = np.loadtxt('ex2data1.txt', delimiter=',')
X = data[:, [0,1]]
y = data[:, 2]

次に,試験1の点数を横軸,試験2の点数を縦軸とした平面上に,合格者についてはx印で,不合格者については印でプロットされた散布図を作成します.これは次のような関数で書けます:

def plot_data(X, y):
    X_1 = X[y == 1] # 合格者に関する試験の点数データ
    X_0 = X[y == 0] # 不合格者に関する試験の点数データ
    plt.scatter(X_1[:, 0], X_1[:, 1], marker='+', label='Admitted')
    plt.scatter(X_0[:, 0], X_0[:,1], marker='o', label='not admitted')
    plt.xlabel('Exam 1 score')
    plt.ylabel('Exam 2 score')
    plt.legend()
    plt.show()

ここでのポイントは,booleanマスクを使って,合格者の試験成績データと不合格者の試験成績データに分けている部分です.

これを用いると次のような散布図を作成できます:
c2.png

損失関数とその微分を計算

今回のロジスティック回帰に関する仮定関数は次で与えられます:

$$
h_{\theta}(x) = g(\theta^Tx),
$$
ここで,$g$はシグモイド関数
$$
g(z) = \frac{1}{1 + \exp(-z)}
$$
です.

まず,シグモイド関数を実装してみましょう.以下が実装例です.

def sigmoid(z):
    g = 1 / (1 + np.exp(-z))
    return g

次に,ロジスティック回帰の場合の損失関数とその微分を思い出しましょう.

損失関数は以下であり,
$$
J(\theta) = \frac{1}{m} \sum_{i=1}^m [ -y^{(i)}\log(h_{\theta}(x^{(i)})) - ( 1 - y^{(i)})\log(1 - h_{\theta}(x^{(i)})) ] \tag{1}
$$
その微分は,
$$
\frac{\partial J(\theta) }{\partial \theta_j} = \frac{1}{m}\sum_{i=1}^m (h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)}\quad (j = 0, 1, \ldots,n) \tag{2}
$$
となります.

式(1)の微分が(2)で与えられることの証明

一旦ここで式(1)の微分が(2)で与えられることの証明を与えてみましょう.直接計算で示せることですので,やるだけです.

まず,式(1)の左辺を$\theta_j$に関して偏微分すると

\begin{align}
\frac{\partial J(\theta)}{\partial \theta_j} &= \frac{1}{m} \sum_{i=1}^m [ -y^{(i)}\frac{\partial}{\partial \theta_j}\log(h_{\theta}(x^{(i)})) - ( 1 - y^{(i)})\frac{\partial}{\partial \theta_j}\log(1 - h_{\theta}(x^{(i)})) ]  \\
&= \frac{1}{m} \sum_{i=1}^m [ -y^{(i)}\frac{1}{h_{\theta}(x^{(i)})}\frac{\partial}{\partial \theta_j}h_{\theta}(x^{(i)}) - ( 1 - y^{(i)})\frac{-1}{1 - h_{\theta}(x^{(i)})}\frac{\partial}{\partial \theta_j}h_{\theta}(x^{(i)}) ] \tag{3}
\end{align}

となります.ここで,

$$
\frac{\partial}{\partial \theta_j}h_{\theta}(x^{(i)}) = \frac{x_j^{(i)}e^{-\theta^Tx}}{(1 + e^{-\theta^T x})^2}, \tag{4}
$$
$$
\frac{1}{h_{\theta}(x^{(i)})} = 1 + e^{- \theta^Tx}, \tag{5}
$$
$$
\frac{1}{1 - h_{\theta}(x^{(i)})} = \frac{1 + e^{-\theta^T x}}{e^{-\theta^Tx}} \tag{6}
$$
となることから,式(3)に式(4),(5),(6)を代入すると,

\begin{align}
\frac{\partial J(\theta)}{\partial \theta_j} &=  \frac{1}{m} \sum_{i=1}^m [ -y^{(i)}(1 + e^{- \theta^Tx})\frac{x_j^{(i)}e^{-\theta^Tx}}{(1 + e^{-\theta^T x})^2} - ( 1 - y^{(i)})\frac{1 + e^{-\theta^T x}}{e^{-\theta^Tx}}\frac{x_j^{(i)}e^{-\theta^Tx}}{(1 + e^{-\theta^T x})^2} ] \\
&= \frac{1}{m}\sum_{i=1}^m[ - y^{(i)}x_j^{(i)}e^{-\theta^T x} + (1 - y^{(i)})\frac{x_j^{(i)}}{1 + e^{-\theta^T x}}]\\
& =\frac{1}{m}\sum_{i=1}^m (h_{\theta}(x^{(i)}) - y^{(i)})x_j^{(i)}
\end{align}

となり,式(2)が示せました.

実装

それでは,損失関数とその微分を実装しましょう.式をコードにすると次のようになります.

def cost_function(theta, *args):
    X, y = args
    h = sigmoid(np.matmul(X, theta))
    J = np.mean(-y * np.log(h + 1e-10) - (1 - y) * np.log((1 - h) + 1e-10))
    return J
def gradient(theta, *args):
    X, y = args
    h = sigmoid(np.matmul(X, theta))
    grad = np.matmul((h - y), X) / m 
    return grad

ここでは,あとでfmin_cgを使う関係で,それに合わせた書き方にしてみました.また$\log(x)$は$x$が小さくなると不安定になるので,1e-10を入れておきました(これがどの程度意味があるのかこの辺りは詳しくないので、詳しい人がいたら教えてほしい).

学習

元々の問題では,Octaveのfminuncを使って学習するようになっていたので,Pythonでやる上では,scipyのfmin_cgを使ってやってみました.

theta = optimize.fmin_cg(cost_function, x0=initial_theta, args=(X,y), fprime=gradient, maxiter=400)

Decision boundaryを描く

この部分のコードは問題にはなっていなかったのですが,Pythonバージョンとして以下のように書き直しました.

def plot_decision_boundary(theta, X, y):
    X_1 = X[y == 1]
    X_0 = X[y == 0]
    plt.scatter(X_1[:, 1], X_1[:,2], marker='+', label='Admitted')
    plt.scatter(X_0[:, 1], X_0[:,2], marker='o', label='not admitted')

    # Only need 2 points to define a line, so choose two endpoints
    plot_x = np.array([np.min(X[:,1]) - 2, np.max(X[:, 1]) + 2])
    # Calculate the decision boundary line
    plot_y = - (theta[0] + theta[1] * plot_x) / theta[2]

    plt.plot(plot_x, plot_y, label='Decision Boundary', color='r')
    plt.axis([30, 100, 30, 100])
    plt.legend()
    plt.show()

これを用いると次のようにdecision boundaryを散布図に追加した図が表示できます:
c2-1.png

予測

応募者の合格確率は0.5を閾値として,仮定関数が0.5以上であれば合格(1),小さければ不合格(0)とします.これはコードでは以下のように書けます:

def predict(theta, X):
    m = X.shape[0]
    h = sigmoid(np.matmul(X, theta))
    p = (h >= 0.5).astype(np.uint8)

    return p

コード全体

今回の問題のコード全体は以下のようになります.

import numpy as np
import matplotlib.pyplot as plt
from scipy import optimize

data = np.loadtxt('ex2data1.txt', delimiter=',')

X = data[:, [0, 1]]
y = data[:, 2]

print('Plotting data with + indicating (y = 1) examples and o indicating (y = 0) examples.')


def plot_data(X, y):
    X_1 = X[y == 1]
    X_0 = X[y == 0]
    plt.scatter(X_1[:, 0], X_1[:, 1], marker='+', label='Admitted')
    plt.scatter(X_0[:, 0], X_0[:, 1], marker='o', label='not admitted')
    plt.xlabel('Exam 1 score')
    plt.ylabel('Exam 2 score')
    plt.legend()
    plt.show()

plot_data(X, y)


def sigmoid(z):
    g = 1 / (1 + np.exp(-z))

    return g

m, n = X.shape
X = np.column_stack((np.ones([m]), X))

initial_theta = np.zeros([n + 1])


def cost_function(theta, *args):
    X, y = args
    h = sigmoid(np.matmul(X, theta))
    J = np.mean(-y * np.log(h + 1e-10) - (1 - y) * np.log((1 - h) + 1e-10))
    return J


def gradient(theta, *args):
    X, y = args
    h = sigmoid(np.matmul(X, theta))
    grad = np.matmul((h - y), X) / m

    return grad

cost, grad = cost_function(initial_theta, *(X, y)), gradient(initial_theta, *(X, y))

print('Cost at initial theta (zeros): {:f}'.format(cost))
print('Expected cost (approx): 0.693')
print('Gradient at initial theta (zeros):')
print(grad)
print('Expected gradients (approx):\n[-0.1000  -12.0092  -11.2628]')

test_theta = np.array([-24, 0.2, 0.2])

cost, grad = cost_function(test_theta, *(X, y)), gradient(test_theta, *(X, y))

print('Cost at test theta: {:f}'.format(cost))
print('Expected cost (approx): 0.218')
print('Gradient at test theta:')
print(grad)
print('Expected gradients (approx):\n[0.043  2.566  2.647')

input('Program paused. Press enter to continue.')

theta = optimize.fmin_cg(cost_function, x0=initial_theta, args=(X, y), fprime=gradient, maxiter=400)

cost = cost_function(theta, *(X, y))

print('Cost at theta found by fminunc: {:f}'.format(cost));
print('Expected cost (approx): 0.203')
print('theta:')
print(theta)
print('Expected theta (approx):')
print('[-25.161  0.206  0.201]')


def plot_decision_boundary(theta, X, y):
    X_1 = X[y == 1]
    X_0 = X[y == 0]
    plt.scatter(X_1[:, 1], X_1[:, 2], marker='+', label='Admitted')
    plt.scatter(X_0[:, 1], X_0[:, 2], marker='o', label='not admitted')

    # Only need 2 points to define a line, so choose two endpoints
    plot_x = np.array([np.min(X[:, 1]) - 2, np.max(X[:, 1]) + 2])
    # Calculate the decision boundary line
    plot_y = - (theta[0] + theta[1] * plot_x) / theta[2]

    plt.plot(plot_x, plot_y, label='Decision Boundary', color='r')
    plt.axis([30, 100, 30, 100])
    plt.legend()
    plt.show()


plot_decision_boundary(theta, X, y)

input('Program paused. Press enter to continue.')

prob = sigmoid(np.matmul([1, 45, 85], theta))

print('For a student with scores 45 and 85, we predict an admission probability of {:f}'.format(prob))
print('Expected value: 0.775 +/- 0.002')


def predict(theta, X):
    m = X.shape[0]
    h = sigmoid(np.matmul(X, theta))
    p = (h >= 0.5).astype(np.uint8)
    return p


p = predict(theta, X)

print('Train Accuracy: {:f}'.format(np.mean(p == y) * 100));
print('Expected accuracy (approx): 89.0')

参考

https://www.coursera.org/learn/machine-learning/