实现了基于线性最小二乘、基于梯度下降法、基于随机梯度下降法的多元线性回归。计算中发现对一元回归，三者均有较高的R方值；对多元回归，基于线性最小二乘的拟合精度最好，后两者精度较差。其原因为梯度下降法由于起始点选择的原因可能陷入局部最优。
源程序
import numpy as np
from sklearn import datasets
from sklearn.metrics import r2_score
class LinearRegression:
def __init__(self):
"""初始化Linear Regression模型"""
self.coef_ = None
self.intercept_ = None
self._theta = None
def scaling(self, X_train, predict_X):
# """最值归一化"""
# scaling_X = np.ones(shape=original_X.shape)
# for i in range(original_X.shape[1]):
# scaling_X[:,i] = (original_X[:,i] - np.min(X_train[:,i]))/(np.max(X_train[:,i])-np.min(X_train[:,i]))
"""均值方差归一化"""
scaling_X = np.empty(shape=X_train.shape)
mean_X = np.empty(shape=X_train.shape[1])
std_X = np.empty(shape=X_train.shape[1])
for i in range(X_train.shape[1]):
scaling_X[:,i] =(predict_X[:,i] - np.mean(X_train[:,i])) / (np.std(X_train[:,i]))
mean_X[i] = np.mean(X_train[:,i])
std_X[i] = np.std(X_train[:,i])
return scaling_X, mean_X, std_X
def fit_normal(self, X_train, y_train):
"""根据最小二乘法训练Linear Regression模型"""
assert X_train.shape[0] == y_train.shape[0], \
"the size of X_train must be equal to the size of y_train"
X_b = np.hstack([np.ones((len(X_train), 1)), X_train]) #自变量左边加上全为1的一列
self._theta = np.linalg.inv(X_b.T.dot(X_b)).dot(X_b.T).dot(y_train) #计算拟合的函数值
self.intercept_ = self._theta[0].copy() #截距
self.coef_ = self._theta[1:].copy() #系数
return self
def fit_gd(self, X_train, y_train, eta=0.01, n_iters=1e4):
"""根据训练数据集X_train, y_train, 使用梯度下降法训练Linear Regression模型"""
assert X_train.shape[0] == y_train.shape[0], \
"the size of X_train must be equal to the size of y_train"
def J(theta, X_b, y):
try:
return np.sum((y - X_b.dot(theta)) ** 2) / len(y)
except:
return float('inf')
def dJ(theta, X_b, y):
# res = np.empty(len(theta))
# res[0] = np.sum(X_b.dot(theta) - y)
# for i in range(1, len(theta)):
# res[i] = (X_b.dot(theta) - y.reshape(-1)).dot(X_b[:, i])
# return res * 2 / len(X_b)
return (X_b.T).dot(X_b.dot(theta) - y.reshape(-1)) * 2 / len(y) #X_b.dot(theta)计算结果为一个一维向量
def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):
theta = initial_theta.copy() #初始化拟合参数向量
cur_iter = 0 #迭代循环次数
while cur_iter < n_iters:
gradient = dJ(theta, X_b, y)
last_theta = theta.copy()
theta = theta - eta * gradient
if (abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):
break
cur_iter += 1
return theta
X_train, mean_X_train, std_X_train= self.scaling(X_train, X_train)
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
initial_theta = np.zeros(X_b.shape[1])
self._theta = gradient_descent(X_b, y_train, initial_theta, eta, n_iters)
#将所拟合的参数做归一化的逆向运算
self.intercept_ = self._theta[0] - np.sum(self._theta[1:] * mean_X_train / std_X_train)
self.coef_ = self._theta[1:] / std_X_train
return self
def fit_sgd(self, X_train, y_train, n_iters=5, t0=5, t1=50):
"""根据训练数据集X_train, y_train, 使用随机梯度下降法训练Linear Regression模型"""
assert X_train.shape[0] == y_train.shape[0], \
"the size of X_train must be equal to the size of y_train"
assert n_iters >= 1
def dJ_sgd(theta, X_b_i, y_i):
return X_b_i * (X_b_i.dot(theta) - y_i) * 2.
def sgd(X_b, y, initial_theta, n_iters, t0=5, t1=50):
def learning_rate(t):
return t0 / (t + t1)
theta = initial_theta.copy()
m = len(X_b)
for cur_iter in range(n_iters):
indexes = np.random.permutation(m)
X_b_new = X_b[indexes].copy()
y_new = y[indexes].copy()
for i in range(m):
gradient = dJ_sgd(theta, X_b_new[i], y_new[i])
theta = theta - learning_rate(cur_iter * m + i) * gradient
return theta
X_train, mean_X_train, std_X_train= self.scaling(X_train, X_train)
X_b = np.hstack([np.ones((len(X_train), 1)), X_train])
initial_theta = np.random.randn(X_b.shape[1])
self._theta = sgd(X_b, y_train, initial_theta, n_iters, t0, t1)
#将所拟合的参数做归一化的逆向运算
self.intercept_ = self._theta[0] - np.sum(self._theta[1:] * mean_X_train / std_X_train)
self.coef_ = self._theta[1:] / std_X_train
return self
def predict(self, X_predict):
"""给定待预测数据集X_predict，返回表示X_predict的结果向量"""
assert self.intercept_ is not None and self.coef_ is not None, \
"must fit before predict!"
assert X_predict.shape[1] == len(self.coef_), \
"the feature number of X_predict must be equal to X_train"
X_b = np.hstack([np.ones((len(X_predict), 1)), X_predict])
return X_b.dot(self._theta)
def score_normal(self, X_test, y_test):
"""根据测试数据集 X_test 和 y_test 确定最小二乘模型的准确度"""
y_predict = self.predict(X_test)
return r2_score(y_test, y_predict)
def score_gd(self, X_train, X_test, y_test):
"""根据测试数据集 X_test 和 y_test 确定梯度下降、随机梯度下降模型的准确度"""
X_test, mean_X_train, std_X_train= self.scaling(X_train, X_test)
y_predict = self.predict(X_test)
return r2_score(y_test, y_predict)
def __repr__(self):
return "LinearRegression()"
#boston = datasets.load_boston()
##print (boston.DESCR)
#X_train = boston.data
#y_train = boston.target
#X_train = X_train[y_train<50]
#y_train = y_train[y_train<50]
X_train = np.random.normal(3,4,size = (100,3))
y_train = X_train.dot(np.array([[1], [2], [3]])) + 3 + np.random.normal(0,1,size = (100,1))
Linearmodel = LinearRegression()
Linearmodel.fit_normal(X_train, y_train)
print(Linearmodel.score_normal(X_train, y_train))
print(Linearmodel.intercept_, Linearmodel.coef_, '\\n')
Linearmodel.fit_gd(X_train, y_train, eta=1e-3, n_iters=1e5)
print(Linearmodel.score_gd(X_train, X_train, y_train))
print(Linearmodel.intercept_, Linearmodel.coef_, '\\n')
Linearmodel.fit_sgd(X_train, y_train, n_iters=500, t0=5, t1=50)
print(Linearmodel.score_gd(X_train, X_train, y_train))
print(Linearmodel.intercept_, Linearmodel.coef_, '\\n')
python tips
X_b.dot(theta) - y.reshape(-1)的计算中X_b.dot(theta)计算结果为一个1001一维向量,y为一个1001的矩阵，虽然只1列，但仍为二维矩阵，故两者相减将报错。
对于使用基于梯度下降法、基于随机梯度下降法的多元线性回归均需要进行数据归一化，而基于最小二乘的线性回归不需要。
np.hstack([A, B])将A、B两矩阵水平方向拼合；np.vstack([A, B])将A、B两矩阵竖直方向拼合；np.concatenate([A, B],axis=1)通过制定axis值控制水平和竖直方向的拼接。
函数shuffle与permutation都是对原来的数组进行重新洗牌（即随机打乱原来的元素顺序）；区别在于shuffle直接在原来的数组上进行操作，改变原来数组的顺序，无返回值。而permutation不直接在原来的数组上进行操作，而是返回一个新的打乱顺序的数组，并不改变原来的数组。
