「Machine Learning」(第2週)の復習① (2017/06/21)

今回は第2週の分の前半である。プログラムの実行は課題の「ex2.m」ファイルの順序に基づいている。

準備

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

# Machine Learning Online Class - Exercise 2: Logistic Regression

# Initialization

# Load Data

data = np.loadtxt("ex2data1.txt",delimiter=",")
x = data[:, 0:2]
y = data[:, 2]

Part 1: Plotting

# ==================== Part 1: Plotting ====================

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

def plot_data(x, y):
    pos = np.where(y==1)
    neg = np.where(y==0)
    plt.plot(x[pos, 0].T, x[pos, 1].T, 'k+', linewidth=2, markersize=7,
             label='Admitted')
    plt.plot(x[neg, 0].T, x[neg, 1].T, 'ko', markerfacecolor='y', markersize=7,
             label='Not admitted')

plot_data(x, y)

plt.xlabel('Exam 1 score')
plt.ylabel('Exam 2 score')

plt.legend(loc='upper right')
plt.savefig('ex2_1.png', dpi=150)  # Added.
plt.show()

print('\nProgram paused. Press enter to continue.')

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

画像

Program paused. Press enter to continue.

Part 2: Compute Cost and Gradient

ロジスティック回帰の仮説の数式は以下のようである。

\[ h_{\theta}(x) = g(\theta^{T}x) \]

  • \( \theta \): パラメータ
  • \( x \): 説明変数

これは以下のcost_function()関数内の

z = np.dot(x, theta)
h = sigmoid(z)

に対応している。

また、シグモイド関数の数式は以下のようである。

\[ g(z) = \frac {1} {1 + e^{-z}} \]

  • \( z \): 上述の\( \theta^{T}x \)

これは以下のsigmoid()関数内の

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

に対応している。

さらに、コスト関数の数式は以下のようである。

\[ 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)}))}] \]

  • \( \theta \): パラメータ
  • \( m \): データ数
  • \( x \): 説明変数(正確にはxとθの内積)
  • \( y \): 目的変数

これは以下のcost_function()関数内の

j = 1 / m * np.sum(-y*np.log(h)-(1-y)*np.log(1-h))

に対応している。

最後にコスト関数を偏微分した数式は以下のようである。

\[ \frac {\partial J(\theta)} {\partial \theta_{j}} = \frac {1} {m} \sum_{i=1}^{m} (h_{\theta}(x^{(i)})-y^{(i)}){x_{j}}^{(i)} \]

これは以下のcost_function()関数内の

grad[0, 0] = 1 / m * np.sum((h-y)*x[:, 0].reshape(len(x[:, 0]), 1))
grad[1, 0] = 1 / m * np.sum((h-y)*x[:, 1].reshape(len(x[:, 1]), 1))
grad[2, 0] = 1 / m * np.sum((h-y)*x[:, 2].reshape(len(x[:, 2]), 1))

に対応している。

# ============ Part 2: Compute Cost and Gradient ============
m, n = x.shape

x = np.c_[np.ones([m, 1]), x]

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

def sigmoid(z):
    g = np.zeros(z.shape)

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

def cost_function(theta, x, y):
    y = y.reshape(len(y), 1)
    m = len(y)

    grad = np.zeros(theta.shape)
    z = np.dot(x, theta)
    h = sigmoid(z)

    j = 1 / m * np.sum(-y*np.log(h)-(1-y)*np.log(1-h))
    grad[0, 0] = 1 / m * np.sum((h-y)*x[:, 0].reshape(len(x[:, 0]), 1))
    grad[1, 0] = 1 / m * np.sum((h-y)*x[:, 1].reshape(len(x[:, 1]), 1))
    grad[2, 0] = 1 / m * np.sum((h-y)*x[:, 2].reshape(len(x[:, 2]), 1))
    return j, grad

cost, grad = cost_function(initial_theta, x, y)

print('Cost at initial theta (zeros): ', cost)
print('Expected cost (approx): 0.693')
print('Gradient at initial theta (zeros): ')
print(' ', grad)
print('Expected gradients (approx):\n -0.1000\n -12.0092\n -11.2628')

test_theta = np.array([[-24], [0.2], [0.2]])
cost, grad = cost_function(test_theta, x, y)

print('\nCost at test theta: ', cost)
print('Expected cost (approx): 0.218')
print('Gradient at test theta: ')
print(' ', grad)
print('Expected gradients (approx):\n 0.043\n 2.566\n 2.647')

print('\nProgram paused. Press enter to continue.')

Cost at initial theta (zeros):  0.69314718056
Expected cost (approx): 0.693
Gradient at initial theta (zeros): 
  [[ -0.1       ]
 [-12.00921659]
 [-11.26284221]]
Expected gradients (approx):
 -0.1000
 -12.0092
 -11.2628

Cost at test theta:  0.218330193827
Expected cost (approx): 0.218
Gradient at test theta: 
  [[ 0.04290299]
 [ 2.56623412]
 [ 2.64679737]]
Expected gradients (approx):
 0.043
 2.566
 2.647

Program paused. Press enter to continue.

Part 3: Optimizing using fminunc

# ============= Part 3: Optimizing using fminunc  =============

#options = optimset('GradObj', 'on', 'MaxIter', 400);  # removed

def cost_function2(theta, x, y):
    theta = theta.reshape(len(theta), 1)
    y = y.reshape(len(y), 1)
    m = len(y)

    grad = np.zeros(theta.shape)
    z = np.dot(x, theta)
    h = sigmoid(z)

    j = 1 / m * np.sum(-y*np.log(h)-(1-y)*np.log(1-h))
    grad[0, 0] = 1 / m * np.sum((h-y)*x[:, 0].reshape(len(x[:, 0]), 1))
    grad[1, 0] = 1 / m * np.sum((h-y)*x[:, 1].reshape(len(x[:, 1]), 1))
    grad[2, 0] = 1 / m * np.sum((h-y)*x[:, 2].reshape(len(x[:, 2]), 1))
    return j

result = minimize(cost_function2, initial_theta, args=(x, y), method='TNC')
cost = result.fun
theta = result.x

print('Cost at theta found by fminunc: ', cost)
print('Expected cost (approx): 0.203')
print('theta: ')
print('  ', theta)
print('Expected theta (approx):')
print(' -25.161\n 0.206\n 0.201')

def map_feature(x1, x2):
    degree = 6
    out = np.ones(x1[:, 0].shape)
    end = len(out)  # Added.
    for i in range(degree):
        for j in range(i):
            out[:, end+1] = [x1**(i-j)*(x2**j)]
    return out

def plot_decision_boundary(theta, x, y):
    plot_data(x[:, 1:3], y)
    if x.shape[1] <= 3:
        plot_x = np.array([np.min(x[:, 1])-2, np.max(x[:, 1])+2])

        plot_y = np.array(-1/theta[2])*(theta[1]*plot_x + theta[0])

        plt.plot(plot_x, plot_y)
        plt.legend()  # Modified.

        plt.xlim([30, 100])
        plt.ylim([30, 100])
    else:
        u = np.linspace(-1, 1.5, 50)
        v = np.linspace(-1, 1.5, 50)

        z = np.zeros(len(u), len(v))
        for i in range(u):
            for j in range(v):
                z[i, j] = map_feature(u[i], v[j])*theta
        z = z.T

        plt.contour(u, v, z, [0, 0], linewidth=2)

plot_decision_boundary(theta, x, y)

plt.xlabel('Exam 1 score')
plt.ylabel('Exam 2 score')

plt.legend(loc='upper right')  # Modified.
plt.savefig('ex2_2.png', dpi=150)  # Added.
plt.show()

print('\nProgram paused. Press enter to continue.')

Cost at theta found by fminunc:  0.203601040162
Expected cost (approx): 0.203
theta: 
   [-26.00633023   0.21305365   0.20823975]
Expected theta (approx):
 -25.161
 0.206
 0.201

画像

Program paused. Press enter to continue.

PythonでOctaveのfminunc()関数に相当するのはSciPyのoptimize.minimize()関数だろう。そこでminimize()関数で代用した。

minimize()関数はリストを戻り値とする関数を引数に取れないようだ。このため、コストとパラメータを返すcost_function()関数を引数に取ろうとするとエラーになる。そこでコストのみを返すcost_function2()関数を作成した。

minimize()関数が返すパラメータはfminunc()関数が返すそれと微妙に異なる。設定をいじれば同じ値を返すようにできるかもしれない。ただ、面倒なのでこのままとする。

Part 4: Predict and Accuracies

# ============== Part 4: Predict and Accuracies ==============

prob = sigmoid(np.dot([1, 45, 85], theta))
print('For a student with scores 45 and 85, we predict an admission '
      + 'probability of ', prob)
print('Expected value: 0.775 +/- 0.002\n')

def predict(theta, x):

    m = x.shape[0]

    p = np.zeros([m, 1])

    z = np.dot(x, theta)
    g = sigmoid(z)
    p = g >= 0.5
    return p

p = predict(theta, x)

print('Train Accuracy: ', np.mean((p == y).astype(float)) * 100)
print('Expected accuracy (approx): 89.0')

For a student with scores 45 and 85, we predict an admission probability of  0.782698660738
Expected value: 0.775 +/- 0.002

Train Accuracy:  89.0
Expected accuracy (approx): 89.0

使用したパラメータが微妙に違うため、今回の結果も微妙に異なる。