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

今回は第2週の分の後半である。プログラムの実行は課題の「ex2_reg.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("ex2data2.txt",delimiter=",")
x = data[:, 0:2]
y = data[:, 2]

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='y = 1')
    plt.plot(x[neg, 0].T, x[neg, 1].T, 'ko', markerfacecolor='y', markersize=7,
             label='y = 0')

plot_data(x, y)

plt.xlabel('Microchip Test 1')
plt.ylabel('Microchip Test 2')

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

画像

「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

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

「Machine Learning」(第1週)の復習 (2017/06/09)

Courseraの「Machine Learning」の課題を復習する。

Machine Learning

ただ、課題のコードはOctaveではなくてPythonで書く。課題のコードを公開することは禁じられているが、Pythonで書くなら問題はないだろう(本当か?)。

また、データファイルは課題で与えられたものを借用した。プログラムを実行するときはデータファイルを予め作業ディレクトリに保存しておく。

データファイルは私のGitHubに置いてある。だが、問題があれば削除する。

https://github.com/fxst24/fxst24/tree/master/machine_learning

課題では関数ごとにファイルが作られるが、ここではファイルは作らず、関数はコマンド内で定義している。

今回は第1週の分である。プログラムの実行は課題の「ex1.m」ファイルの順序に基づいている。また、オプションの課題はやっていないので後日追加する。

準備

import numpy as np
import matplotlib.pyplot as plt

# Machine Learning Online Class - Exercise 1: Linear Regression

# Initialization

Part 1: Basic Function

# ==================== Part 1: Basic Function ====================
print('Running warmUpExercise ... ')
print('5x5 Identity Matrix: ')

def warm_up_exercise():
    a = np.identity(5)
    return a
print(warm_up_exercise())

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

Running warmUpExercise ... 
5x5 Identity Matrix: 
[[ 1.  0.  0.  0.  0.]
 [ 0.  1.  0.  0.  0.]
 [ 0.  0.  1.  0.  0.]
 [ 0.  0.  0.  1.  0.]
 [ 0.  0.  0.  0.  1.]]
Program paused. Press enter to continue.

Part 2: Plotting

# ======================= Part 2: Plotting =======================
print('Plotting Data ...')
data = np.loadtxt("ex1data1.txt",delimiter=",")
x = data[:, 0]
y = data[:, 1]
m = len(y)

def plot_data(x, y):
    plt.plot(x, y, 'rx', markersize=10)
    plt.ylabel('Profit in $10,000s')
    plt.xlabel('Population of City in 10,000s')
    plt.savefig('ex1_1.png', dpi=150)  # Added.
    plt.show()

plot_data(x, y)

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

Plotting Data ...

画像

Program paused. Press enter to continue.

Part 3: Cost and Gradient descent

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

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

  • \( \theta \): パラメータ
  • \( m \): データ数
  • \( x \): 説明変数
  • \( y \): 目的変数

これは以下のcompute_cost()関数に対応している。

また、線形モデルの数式は以下のようである。

\[ h_{\theta}(x) = \theta^{T}x = \theta_{0} + \theta_{1}x_{1} \]

これはcompute_cost()関数内の「np.dot(x, theta)」に対応している。

数式ではθの転置行列とxとの内積であるのに、コマンドではxとθとの内積になっているのはおかしいと思うだろう。

データはxが97行2列、θが2行1列になっている。したがって、内積によって97行1列となり、実際これで構わない(この後、総和を求めるので1行97列でも構わない)。

数式のように計算すると、シータの転置行列が1行2列、xが97行2列であるから、これでは形が合わず、内積はできない。数式のように計算できるためにはxが2行97列でなければならないのである。そうすれば1行97列になる。

だが、与えられたデータの形式がこうなっているので、細かいことは気にせず(おいっ!)、臨機応変に対応する。これ以降の課題も同様である。

# =================== Part 3: Cost and Gradient descent ===================
x = np.c_[np.ones([m, 1]), data[:, 0]]
theta = np.zeros([2, 1])

iterations = 1500
alpha = 0.01

print('\nTesting the cost function ...')

def compute_cost(x, y, theta):
    m = len(y)
    y = y.reshape(len(y), 1)
    j = 1 / (2*m) * np.sum((np.dot(x, theta)-y)**2)
    return j

j = compute_cost(x, y, theta)
print('With theta = [0 ; 0]\nCost computed = ', j)
print('Expected cost value (approx) 32.07')

j = compute_cost(x, y, [[-1], [2]])
print('\nWith theta = [-1 ; 2]\nCost computed = ', j)
print('Expected cost value (approx) 54.24')

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

Testing the cost function ...
With theta = [0 ; 0]
Cost computed =  32.0727338775
Expected cost value (approx) 32.07

With theta = [-1 ; 2]
Cost computed =  54.242455082
Expected cost value (approx) 54.24
Program paused. Press enter to continue.

NumPyの1次元配列は行ベクトルでも列ベクトルでもないようだ。だが、行列演算に使うと行ベクトルとして扱われているっぽい。このため、注意が必要だ。

例えば、

y.shape

とすれば、

(97,)

と返ってくる。これは「(97, 1)」でも「(1, 97)」でもないが、行列演算では「(1, 97)」として扱われるようだ。「(97, 1)」として扱いたい場合は

y = y.reshape(len(y), 1)
y.shape

(97, 1)

のようにしてリシェイプする。

さて、勾配降下法の数式は以下のようである。

\[ \theta_{j} := \theta_{j} - \alpha \frac {1} {m} \sum_{i=1}^{m} (h_{\theta}(x^{(i)}) - y^{(i)})x^{(i)}_{j} \]

  • \( \alpha \): 学習率

これは以下のgradient_descent()関数に対応している。

print('\nRunning Gradient Descent ...')

def gradient_descent(x, y, theta, alpha, num_iters):
    m = len(y)
    y = y.reshape(len(y), 1)
    j_history = np.zeros([num_iters, 1])
    for i in range(num_iters):
        temp = theta[0, 0] - alpha * 1 / m * np.sum(np.dot(x, theta)-y)
        theta[1, 0] -= (
                alpha * 1 / m * np.sum((np.dot(x, theta)-y)
                * x[:, 1].reshape(len(x[:, 1]), 1)))
        theta[0, 0] = temp
        j_history[i, 0] = compute_cost(x, y, theta)
    return theta

theta = gradient_descent(x, y, theta, alpha, iterations)

print('Theta found by gradient descent:')
print(theta)
print('Expected theta values (approx)')
print(' -3.6303\n  1.1664\n')

plt.plot(x[:, 1], y, 'rx', markersize=10, label='Training data')  # Added.
plt.plot(x[:, 1], np.dot(x, theta), '-', label='Linear regression')
plt.legend()
plt.savefig('ex1_2.png', dpi=150)  # Added.
plt.show()

predict1 = np.dot([1, 3.5], theta)
print('For population = 35,000, we predict a profit of ', predict1*10000)
predict2 = np.dot([1, 7], theta)
print('For population = 70,000, we predict a profit of ', predict2*10000)

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

Running Gradient Descent ...
Theta found by gradient descent:
[[-3.63029144]
 [ 1.16636235]]
Expected theta values (approx)
 -3.6303
  1.1664

画像

For population = 35,000, we predict a profit of  [ 4519.7678677]
For population = 70,000, we predict a profit of  [ 45342.45012945]
Program paused. Press enter to continue.

xのような2次元配列の1列を切り取った場合、当然、列ベクトルだと思うだろう。ところがそうではない。例えば、

x.shape

(97, 2)

だが、

x[:, 1].shape

(97,)

なのだ。だから、列ベクトルとして扱いたい場合は面倒だが、

x_col1 = x[:, 1].reshape(len(x[:, 1]), 1)
x_col1.shape

(97, 1)

のようにする。しかし、よく分からん仕様だな。