「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)

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