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

ex1_m.py

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

# おまけ1（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()

# データの標準化
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)

# おまけ2（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()


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

# 収束グラフのプロット
plt.plot(np.arange(len(J_history)), J_history, color="b")
plt.xlabel("Number of iterations")
plt.ylabel("Cost J")
plt.show()

# 最急降下法の結果の表示（スケール調整しているので調整前とは結果が違うので注意）
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()

# 正規方程式
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)