Coursera Machine LearningをPythonで実装 - [Week2]単回帰分析、重回帰分析

機械学習の学習素材として大人気の、スタンフォード大学のAndrew Ng先生による[Courseraの講義](https://www.coursera.org/learn/machine-learning/)を、復習がてらにPythonで書き換えてみました。講義ではOctave/Matlabを使っていますが、機械学習やるからにはPythonでやりたいですよね。[Honor Code](https://www.coursera.org/about/terms/honorcode)に引っかからない程度で（これをそのままコピペしても課題はできません）。2018年4月時点のものです。

**追記**：gistからソースファイルをダウンロードできます
単回帰分析：https://gist.github.com/koshian2/c26ee4d139fa1083cabeb8aa2a53696b
重回帰分析：https://gist.github.com/koshian2/8b45ff6cf6d65bf30da7935689790402

#単回帰分析
##（おすすめ1）sklearn.linear_modelのLinearRegressionを使う
```python
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

# データ
X = np.array([6.1101,5.5277,8.5186,7.0032,5.8598,8.3829,7.4764,8.5781,6.4862,5.0546,5.7107,14.164,5.734,8.4084,5.6407,5.3794,6.3654,5.1301,6.4296,7.0708,6.1891,20.27,5.4901,6.3261,5.5649,18.945,12.828,10.957,13.176,22.203,5.2524,6.5894,9.2482,5.8918,8.2111,7.9334,8.0959,5.6063,12.836,6.3534,5.4069,6.8825,11.708,5.7737,7.8247,7.0931,5.0702,5.8014,11.7,5.5416,7.5402,5.3077,7.4239,7.6031,6.3328,6.3589,6.2742,5.6397,9.3102,9.4536,8.8254,5.1793,21.279,14.908,18.959,7.2182,8.2951,10.236,5.4994,20.341,10.136,7.3345,6.0062,7.2259,5.0269,6.5479,7.5386,5.0365,10.274,5.1077,5.7292,5.1884,6.3557,9.7687,6.5159,8.5172,9.1802,6.002,5.5204,5.0594,5.7077,7.6366,5.8707,5.3054,8.2934,13.394,5.4369])
y = np.array([17.592,9.1302,13.662,11.854,6.8233,11.886,4.3483,12,6.5987,3.8166,3.2522,15.505,3.1551,7.2258,0.71618,3.5129,5.3048,0.56077,3.6518,5.3893,3.1386,21.767,4.263,5.1875,3.0825,22.638,13.501,7.0467,14.692,24.147,-1.22,5.9966,12.134,1.8495,6.5426,4.5623,4.1164,3.3928,10.117,5.4974,0.55657,3.9115,5.3854,2.4406,6.7318,1.0463,5.1337,1.844,8.0043,1.0179,6.7504,1.8396,4.2885,4.9981,1.4233,-1.4211,2.4756,4.6042,3.9624,5.4141,5.1694,-0.74279,17.929,12.054,17.054,4.8852,5.7442,7.7754,1.0173,20.992,6.6799,4.0259,1.2784,3.3411,-2.6807,0.29678,3.8845,5.7014,6.7526,2.0576,0.47953,0.20421,0.67861,7.5435,5.3436,4.2415,6.7981,0.92695,0.152,2.8214,1.8451,4.2959,7.2029,1.9869,0.14454,9.0551,0.61705])

# Scikit-learnの組み込み最小二乗法
from sklearn.linear_model import LinearRegression
print("LinearRegression on sklearn.linear_model ...")
regr = LinearRegression(normalize=False) #Trueにした場合は平均とL2ノルムでスケール調整が行われる
regr.fit(X.reshape(-1, 1), y) #データが1変数の場合はreshapeするのが大事
print("Intercept:", regr.intercept_) #切片（定数項）
print("Coefficients: ", regr.coef_) #xの係数
predicts = regr.predict(np.array([3.5, 7]).reshape(-1, 1))
print("For population = [35k, 70k], we predict a profit of \n", predicts * 10000)
print()
```
忙しい人向け。これで全部結果が出ます。```SGDRegressor```が勾配降下法を使うのに対して、```LinearRegression```は正規方程式で計算します。
正規方程式は一発で解が出て便利ですが、逆行列の計算にO(N^3)のオーダーがかかるので、万、十万単位のデータを扱うときなど、データ数が多いときは最急降下法を使ったほうが速くなる、とのことです（講義より）。

```
LinearRegression on sklearn.linear_model ...
Intercept: -3.89578087831
Coefficients:  [ 1.19303364]
For population = [35k, 70k], we predict a profit of
 [  2798.36876352  44554.54631015]
```

##（おすすめ2）sklearn.linear_modelのSGDRegressorを使う
```python
# Scikit-learnの組み込み勾配降下法による回帰
from sklearn.linear_model import SGDRegressor
print("SGDRegressor on sklearn.linear_model ...")
#shuffle=Falseにすると最急降下法になるはず。max_iterはScikit-learn0.19以降
regr = SGDRegressor(eta0=0.01, learning_rate="constant", max_iter=1500, shuffle=False)
regr.fit(X.reshape(-1, 1), y)
print("Intercept:", regr.intercept_)
print("Coefficients: ", regr.coef_)
predicts = regr.predict(np.array([3.5, 7]).reshape(-1, 1))
print("For population = [35k, 70k], we predict a profit of \n", predicts * 10000)
print()
```

この```SGDRegressor```は確率的勾配降下法なのですが、シャッフルをしなければ最急降下法に近い処理になるので、最急降下法としても使えます。自分で実装する場合と設定値は同じになるよう努めましたが、内部の実装が違うため答えの値は変わります。設定値を変えていろいろ試してみてください。

```
SGDRegressor on sklearn.linear_model ...
Intercept: [-3.84811446]
Coefficients:  [ 1.05704926]
For population = [35k, 70k], we predict a profit of
 [ -1484.42048023  35512.30359183]
```

##自分で実装する
暇な人向け。自分で確認したい場合の方法。1回やれば十分なので、あとは上記の組み込みの方法に任せましょう。

###1.データをプロットする
```python
# プロット
plt.plot(X, y, "o", color="r")
plt.show()
```

![ml_wk2_1.png](https://qiita-image-store.s3.amazonaws.com/0/232088/88c06402-02f4-78e4-d0a4-e35af795168c.png)

###2.コスト関数
```python
# コスト関数
def compute_cost(X, y, theta):
    return np.sum((np.dot(X, theta) - y) ** 2) / 2 / len(y)

print("Testing the cost function ...")
X = np.array((np.ones(len(y)), X)).T
theta = np.zeros(2)
J = compute_cost(X, y, theta)
print("With theta = [0, 0]\nCost computed =", J)
J = compute_cost(X, y, np.array([-1,2]))
print("With theta = [-1, 2]\nCost computed =", J)
print()
```

```
Testing the cost function ...
With theta = [0, 0]
Cost computed = 32.0727338775
With theta = [-1, 2]
Cost computed = 54.242455082
```

###3.最急降下法
```python
def gradient_descent(X, y, theta, alpha, num_iters):
    J_history = []
    theta_out = theta
    for i in range(num_iters):
        theta_out -= alpha / len(y) * np.dot(X.T, np.dot(X, theta_out) - y)
        J_history.append(compute_cost(X, y, theta_out))
    return theta_out, J_history

print("Running Gradient Descent ...")
iterations = 1500;
alpha = 0.01;
theta, J_history = gradient_descent(X, y, theta, alpha, iterations)
print("Theta found by gradient descent:", theta)
```

```
Running Gradient Descent ...
Theta found by gradient descent: [-3.63029144  1.16636235]
```

###4.近似直線のプロット
```python
# 近似直線のプロット
plt.plot(X[:,1], y, "o", label="Training data", color="r")
plt.plot(X[:,1], np.dot(X, theta), label="Linear regression", color="b")
plt.legend()
plt.show()
```
![ml_wk2_2.png](https://qiita-image-store.s3.amazonaws.com/0/232088/d4164a71-5f48-7beb-95c5-fdac6eea2eae.png)

###5.推定値
```python
# xが 35,000 70,000のときの推定値
predict1 = np.dot(np.array([1, 3.5]), theta)
print("For population = 35,000, we predict a profit of ", predict1 * 10000)
predict2 = np.dot(np.array([1, 7]), theta)
print("For population = 70,000, we predict a profit of ", predict2 * 10000)
```

```
For population = 35,000, we predict a profit of  4519.7678677
For population = 70,000, we predict a profit of  45342.4501294
```

###6.誤差関数をプロット
```python
# Jのグラフ化
theta0_vals = np.linspace(-10, 10, 100)
theta1_vals = np.linspace(-1, 4, 100)
theta0, theta1 = np.meshgrid(theta0_vals, theta1_vals)
J_vals = np.zeros((len(theta0_vals), len(theta1_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)
# Surface Plot
fig = plt.figure()
ax = Axes3D(fig)
ax.plot_surface(theta0, theta1, J_vals, rstride=1, cstride=1, cmap=cm.coolwarm)
plt.show()
print()
```

![ml_wk2_3.png](https://qiita-image-store.s3.amazonaws.com/0/232088/2a05d1df-e0fa-c29f-3697-fb7e013a15eb.png)

#重回帰分析
同様に重回帰分析も行います。
##（おすすめ1）sklearn.linear_modelのLinearRegressionを使う
```python
import numpy as np
import matplotlib.pyplot as plt

# データ
X_data = np.array([[2104,3],[1600,3],[2400,3],[1416,2],[3000,4],[1985,4],[1534,3],[1427,3],[1380,3],[1494,3],[1940,4],[2000,3],[1890,3],[4478,5],[1268,3],[2300,4],[1320,2],[1236,3],[2609,4],[3031,4],[1767,3],[1888,2],[1604,3],[1962,4],[3890,3],[1100,3],[1458,3],[2526,3],[2200,3],[2637,3],[1839,2],[1000,1],[2040,4],[3137,3],[1811,4],[1437,3],[1239,3],[2132,4],[4215,4],[2162,4],[1664,2],[2238,3],[2567,4],[1200,3],[852,2],[1852,4],[1203,3]])
y = np.array([399900,329900,369000,232000,539900,299900,314900,198999,212000,242500,239999,347000,329999,699900,259900,449900,299900,199900,499998,599000,252900,255000,242900,259900,573900,249900,464500,469000,475000,299900,349900,169900,314900,579900,285900,249900,229900,345000,549000,287000,368500,329900,314000,299000,179900,299900,239500])

# Scikit-learnの組み込み最小二乗法
from sklearn.linear_model import LinearRegression
print("LinearRegression on sklearn.linear_model ...")
regr = LinearRegression(normalize=False)
regr.fit(X_data, y)
print("Intercept:", regr.intercept_) #切片（定数項）
print("Coefficients: ", regr.coef_) #xの係数
print("Predicted price of a 1650 sq-ft, 3 br house  ... (using LinearRegression):")
predict = regr.predict(np.array([[1650, 3]]))
print("$", predict)
print()
```

使い方は単回帰のときと全く変わりません。

```
LinearRegression on sklearn.linear_model ...
Intercept: 89597.9095428
Coefficients:  [  139.21067402 -8738.01911233]
Predicted price of a 1650 sq-ft, 3 br house  ... (using LinearRegression):
$ [ 293081.4643349]
```

##（おすすめ2）sklearn.linear_modelのSGDRegressorを使う
```python
# データの標準化
def feature_normalize(X):
    mu = np.mean(X, axis=0)
    sigma = np.std(X, axis=0)
    X_norm = (X-mu)/sigma
    return X_norm, mu, sigma

X, mu, sigma = feature_normalize(X_data)

# Scikit-learnの組み込み勾配降下法による回帰
from sklearn.linear_model import SGDRegressor
print("SGDRegressor on sklearn.linear_model ...")
regr = SGDRegressor(eta0=0.01, learning_rate="constant", max_iter=400, shuffle=False)
regr.fit(X, y)
print("Intercept:", regr.intercept_)
print("Coefficients: ", regr.coef_)
predict = regr.predict((np.array([[1650, 3]] - mu) / sigma))
print("Predicted price of a 1650 sq-ft, 3 br house  ... (using SGDRegressor):")
print("$", predict)
print()
```
最急降下法で求めるときは、データのスケールをカラム単位で調整すると収束が速くなるそうです。これはお任せでやる場合も、自分で実装する場合もこれは変わりません。スケール調整する場合は、推定値を求めるときに、代入するパラメーターのスケールを調整してから代入しないといけないので注意しましょう。

```
SGDRegressor on sklearn.linear_model ...
Intercept: [ 339319.80384871]
Coefficients:  [ 106385.95873423   -7991.38807423]
Predicted price of a 1650 sq-ft, 3 br house  ... (using SGDRegressor):
$ [ 293673.79997156]
```

スケール調整済みなので、スケールを調整せずに```LinearRegression```で求めた場合の係数とこの係数では、意味合いが違います。具体的な値を代入して推定値を求めると、```LinearRegression```の場合と似たような結果が得られることを確認できます。

##自分で実装する
データのスケール調整は既に行っているものとします。```SGDRegressor```の例の続きからの想定です。

###1.最急降下法
```python
# Xに切片を追加
X = np.c_[np.ones(len(X)), X]

# コスト関数
def compute_cost_multi(X, y, theta):
    return np.sum((np.dot(X, theta) - y) ** 2) / 2 / len(y)

# 最急降下法
def gradient_descent_multi(X, y, theta, alpha, num_iters):
    J_history = []
    theta_out = theta
    for i in range(num_iters):
        theta_out -= alpha / len(y) * np.dot(X.T, np.dot(X, theta_out) - y)
        J_history.append(compute_cost_multi(X, y, theta_out))
    return theta_out, J_history

print("Running gradient descent ...")

alpha = 0.01
num_iters = 400
theta = np.zeros(3)
theta, J_history = gradient_descent_multi(X, y, theta, alpha, num_iters)
```

###2.収束のプロット
```python
# 収束グラフのプロット
plt.plot(np.arange(len(J_history)), J_history, color="b")
plt.xlabel("Number of iterations")
plt.ylabel("Cost J")
plt.show()
```
![ml_wk2_4.png](https://qiita-image-store.s3.amazonaws.com/0/232088/66f4b4db-5707-b65f-3c1a-d406b855fff0.png)

###3.最急降下法の結果
```python
# 最急降下法の結果の表示（スケール調整しているので調整前とは結果が違うので注意）
print("Theta computed from gradient descent:")
print(theta)
print("Predicted price of a 1650 sq-ft, 3 br house  ... (using gradient descent):")
params = np.append(1, (np.array([1650, 3]) - mu) / sigma)
price = np.dot(theta, params.T)
print("$", price)
print()
```

```
Running gradient descent ...
Theta computed from gradient descent:
[ 334302.06399328   99411.44947359    3267.01285407]
Predicted price of a 1650 sq-ft, 3 br house  ... (using gradient descent):
$ 289221.547371
```

実装が悪かったのか、```SGDRegressor```のときとだいぶ異なる係数が得られました。

###4.正規方程式
```python
# 正規方程式
def normal_eqn(X, y):
    return np.dot(np.dot(np.linalg.pinv(np.dot(X.T, X)), X.T), y)
X = np.c_[np.ones(len(X_data)), X_data]
theta = normal_eqn(X, y)
print("Theta computed from the normal equations:")
print(theta)
print("Predicted price of a 1650 sq-ft, 3 br house  ... (using normal equations):")
price = np.dot(np.array([1, 1650, 3]), theta)
print("$", price)
```

やっていることは```LinearRegression```と同じなので、同一の結果が得られたことに注目（誤差の関係で多少ずれる）。

```
Theta computed from the normal equations:
[ 89597.90954361    139.21067402  -8738.01911255]
Predicted price of a 1650 sq-ft, 3 br house  ... (using normal equations):
$ 293081.464335
```


**次回に進む**
[Coursera Machine LearningをPythonで実装 - [Week3]ロジスティック回帰](https://qiita.com/koshian2/items/eeef29f28dbde1803e6c)