Python
numpy
MachineLearning
matplotlib

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

Courseraの機械学習のコースの宿題はOctaveとMatlabでやるようになっていますが,個人的には普段Octaveは全然使わないので,自分の慣れているPythonでやってみました(TensorFlowとかは使わずNumPyのみでトライ).

Prerequisite

  • Python 3.6.4
  • NumPy 1.13.3
  • Matplotlib 2.0.2

Programming Exercise 1: Linear Regression

Courseraの機械学習のコースの最初の宿題は線形回帰問題です.

今回の問題は,あなたがレストランフランチャイズのCEOであるとして,次にどの都市にお店を出せば良いかを線形回帰で決めようというものです.データとしては都市の人口$X$と,その都市で既にに営業しているレストランでの利益$y$が与えれています.線形回帰で,候補地の人口からそこでの利益を予測します.

宿題を解く

ここでは,ポイントごとに問題を解いたコードを解説します.

最初に以下のように必要なライブラリ・モジュールをインポート:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D

まずデータを読み込もう

データのファイルはex1.txtです.このファイル内の一列目が人口$X$,二列目が売り上げ$y$です.まずはこれをNumPyで読み込んで,人口データからなるベクトル$X$と利益データからなるベクトル$y$に分けます.

NumPyで.txtファイルを読み込むにはnumpy.loadtxtを使います:

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

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

次にこのデータの散布図をmatplotlibで描画します:

plt.scatter(X, y)
plt.show()

すると次のような散布図が出力されます:

ex1.png

損失関数を計算

まず,ベクトル表記で計算するために,人口データに便宜的に1を追加します:

m = len(y)
one = np.ones([m])
X = np.stack((one, X), axis=-1)

これにより,人口の個々のデータは[1.0, 人口]の形になりました.従って,モデルの仮定関数はベクトルを用いた形で以下のように書けます:

$$
h_{\theta}(x) = \theta \cdot x = \theta_0 + \theta_1 x_1
$$

線形回帰問題での損失関数は以下の平均二乗誤差です(ここでは微分した後に簡単になるため1/2を掛けています):

$$
J(\theta) = \frac{1}{2m}\sum_{i=1}^{m} (h_{\theta}(x^{(i)}) - y^{(i)})^2
$$
ここで,$m$はデータ数です.

上の式を素直にコードに落とすと次のようになります:

def compute_cost(X, y, theta):
    # 線形回帰の損失関数を計算する
    # データXとyにフィットする線形回帰のパラメータとしてthetaを用いて損失関数を計算

    h = np.matmul(X, theta)
    J = 0.5 * np.mean((h - y) ** 2)

    return J

勾配降下法のコードを書く

機械学習では,損失関数が最小になるようなパラメータ$\theta$を求めることが重要です.今回の線形回帰では,平均二乗誤差を最小にするような$\theta$がデータにフィットする最適な線形関数を与えます.

今回は,$\theta$を学習させるためのアルゴリズムとして勾配降下法(gradient descent)を用います.これは,以下のように,損失関数の微分を用いて学習率$\alpha$の歩幅で少しパラメータ$\theta$を更新し,それを繰り返すことで損失関数の極小点にたどり着こうというアルゴリズムです:

\begin{align}
\theta_j &:= \theta_j - \alpha \frac{\partial}{\partial \theta_j}J(\theta)\\
&= \theta_j − \alpha \cdot \frac{1}{m}\sum_{i=1}^m (h_{\theta}(x^{(i)}) − y^{(i)})x^{(i)}_j\quad ( j = 0, 1)
\end{align}

ここでの注意点はパラメータの更新は全ての成分を同時に更新する必要があることです.

これをコードで書くと次のようになります:

def gradient_descent(X, y, theta, alpha, num_iters):
    #  thetaを学習させるために勾配降下法を行う
    m = len(y)
    J_history = np.zeros([num_iters])
    for i in range(num_iters):
        h = np.matmul(X, theta)
        grad = np.matmul((h - y), X) / m
        print('iter: {}   gradient : '.format(i))
        print(grad)
        theta = theta - alpha * grad

        cost = compute_cost(X, y, theta)
        print('iter: {}   cost function: {:f}'.format(i, cost))
        #   各イテレーションごとに損失関数Jを保存
        J_history[i] = cost

    return J_history, theta

ここでは,微分の値と損失関数の推移を標準出力に表示するようにし,さらに損失関数のhistoryをリターンするようにしてみました.

学習後のモデルの仮定関数は次のようになります:
ex1_2.png

コード全体

残りはMatplotでグラフを表示する問題なので,コード全体を載せることで解説に代えようと思います.以下が今回の問題のコードとなります:

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import Axes3D


def compute_cost(X, y, theta):

    h = np.matmul(X, theta)
    J = 0.5 * np.mean((h - y) ** 2)

    return J


def gradient_descent(X, y, theta, alpha, num_iters):

    m = len(y)
    J_history = np.zeros([num_iters])
    for i in range(num_iters):
        h = np.matmul(X, theta)
        grad = np.matmul((h - y), X) / m
        print('iter: {}   gradient : '.format(i))
        print(grad)
        theta = theta - alpha * grad

        cost = compute_cost(X, y, theta)
        print('iter: {}   cost function: {:f}'.format(i, cost))

        J_history[i] = cost

    return J_history, theta


print('PlottingData ...')

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

X = data[:, 0]
y = data[:, 1]
m = len(y)

plt.scatter(X, y)
plt.show()

one = np.ones([m])
X = np.stack((one, X), axis=-1)

theta = np.zeros([2])

iterations = 1500
alpha = 0.01

print('Testing the cost function ...')

J = compute_cost(X, y, [-1, 2])
print('With theta = [-1, 2]\nCost computed = {:f}'.format(J))
print('Expected cost value (approx) 54.24')
input('Program paused. Press enter to continue.')
print('Running Gradient Descent ...')
# run gradinet descent
J_history, theta = gradient_descent(X, y, theta, alpha, iterations)

# print theta to screen
print('Theta found by gradient descent:')
print(theta)
print('Expencted theta values (approx)')
print('[-3.6303  1.1664]')

# Plot the linear fit
plt.scatter(X[:, 1], y, label='Training data', color='#1f77b4')
plt.plot(X[:, 1], np.matmul(X, theta), label='Linear regression', color='#ff7f0e')
plt.legend()
plt.show()

# Predict values for population sizes of 35,000 and 70,000
predict1 = np.matmul([1, 3.5], theta)
print('For population = 35,000, we predict a profit of {:f}'.format(predict1 * 10000))

predict2 = np.matmul([1, 7], theta)
print('For population = 70,000, we predict a profit of {:f}'.format(predict2 * 10000))

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

print('Visualizing J(theta_0, theta_1)')

# Grid over which we will calculate J
theta0_vals = np.linspace(-10, 10, 100)
theta1_vals = np.linspace(-1, 4, 100)

# initialize J_vals to matrix of 0's
J_vals = np.zeros([len(theta0_vals), len(theta1_vals)])
# Fill out J_vals
for i in range(len(theta0_vals)):
    for j in range(len(theta1_vals)):
        t = np.array([theta0_vals[i], theta1_vals[j]])
        J_vals[i, j] = compute_cost(X, y, t)

J_vals = np.transpose(J_vals)

fig = plt.figure(figsize=(14, 6))
ax = fig.add_subplot(1, 2, 1, projection='3d')
ax.plot_surface(theta0_vals, theta1_vals, J_vals)
ax.set_xlabel('theta_0')
ax.set_ylabel('theta_1')
plt.show()
plt.close()

# Contour plot
# Plot J_vals as 15 contours spaced logarithmically between 0.01 and 100

plt.contour(theta0_vals, theta1_vals, J_vals, levels=np.logspace(-2, 3, 20))
plt.xlabel('theta_0')
plt.ylabel('theta_1')
plt.plot(theta[0], theta[1], color='r', marker='x', markersize=10, linewidth=2);
plt.show()
plt.close()

参考

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