1.线性回归闭合形式参数求解的原理
如果定义X为m*(n+1)的矩阵,Y为m1的矩阵,θ为(n+1)1维的矩阵,那么在之前的定义中就可以表示为h(x)=Xθ。则代价函数可以表示为J(θ)=1/2(Xθ-y)Т(Xθ-y),J(θ)为凹函数,我们要让其值最小化,只需对该函数求导,然后令导数为0即可求得θ。对其求导后得到XTXθ-XTy,令其等于0,得到θ=(XTX)^-1XT*y。
2.线性回归梯度下降参数求解的原理
我们构造了拟合函数h(θ),并且得到了损失函数J(θ),我们要求得使J(θ)取得最小值的θ,其原理还是求偏导然后使导数为0,我们对J(θ)求导得到(hθ(x) − y) xj,然后可以得到对θj的更新公式
由于数据量较大,所以采用了随机梯度下降,但是准确度相较于批量梯度下降来说会下降。在我的程序里,由于数据采用矩阵形式存储,所以更新过程可以替换为
其中θ为(n+1)1维,X为m(n+1)维,Y为m*1维。梯度为0用损失函数差值小于1e-18来表示,说明这个点是损失函数的极小值点,但并不一定是最小值点。
3.相关文件
4.程序清单
相关包:
from sklearn.preprocessing import StandardScaler from sklearn.model_selection import train_test_split import pandas as pd import numpy as np import os import time from numpy import median from sklearn.preprocessing import OneHotEncoder
(一)读取数据:
# 读取数据 HOUSE_PATH = './' def load_housing_data(housing_path=HOUSE_PATH): csv_path = os.path.join(housing_path, 'housing.csv') return pd.read_csv(csv_path) housing = load_housing_data()
(二)数据处理:
# 将中位数补全空位 median = housing["total_bedrooms"].median() housing["total_bedrooms"].fillna(median, inplace=True) # 独热编码 housing_category = housing[["ocean_proximity"]] cat_encoder = OneHotEncoder() housing_category_onehot = cat_encoder.fit_transform(housing_category) housing = housing.drop("ocean_proximity", axis=1) housing_values = np.c_[housing.values, housing_category_onehot.toarray()] housing_fixed = pd.DataFrame( housing_values, columns=list(housing.columns) + ['<1H OCEAN', 'INLAND', 'ISLAND', 'NEAR BAY', 'NEAR OCEAN'], index=housing.index )
(三)分析数据相关性
# 分析数据相关性 corr_matrix = housing_fixed.corr() # 用corr计算两两特征之间的相关性系数 correlation = corr_matrix["median_house_value"].sort_values(ascending=False) # 跟街区价格中位数特征的其他特征的相关系数 print(correlation)
(四)将数据集分类
# 将数据集分类 train_set, test_set = train_test_split( housing_fixed, test_size=0.3, random_state=42) X_train = train_set[['longitude', 'latitude', 'housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_income', '<1H OCEAN', 'INLAND', 'ISLAND', 'NEAR BAY', 'NEAR OCEAN']] y_train = train_set[["median_house_value"]] X_test = test_set[['longitude', 'latitude', 'housing_median_age', 'total_rooms', 'total_bedrooms', 'population', 'households', 'median_income', '<1H OCEAN', 'INLAND', 'ISLAND', 'NEAR BAY', 'NEAR OCEAN']] y_test = test_set[["median_house_value"]] X = np.hstack([np.ones((len(X_train_std), 1)), X_train_std]) # 训练集X Y = np.array(y_train_std) # 训练集Y x = np.hstack([np.ones((len(X_test_std), 1)), X_test_std]) # 测试集x y = np.array(y_test_std) # 测试集y y_var = np.var(y) # 标准差
(五)特征标准化
# 特征标准化 stdsc = StandardScaler() X_train_std = stdsc.fit_transform(X_train) X_test_std = stdsc.transform(X_test) y_train_std = stdsc.fit_transform(y_train) y_test_std = stdsc.fit_transform(y_test)
(六)初始化θ
theta = np.zeros((14, 1)) # 初始化theta
(七)正规方程
# 正规方程 def nomal(X, Y): theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y) return theta
# 损失函数 def cost_function(x, theta, y): cost = np.sum((np.dot(x, theta)-y)**2) return cost/(2*len(y)) # 正规方程 def nomal(X, Y): theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(Y) return theta # 梯度下降方向 def gradient(X, theta, Y): return X.T.dot((X.dot(theta)-Y))/len(Y) # 梯度下降 def gradient_descent(X, theta, Y, eta): while True: last_theta = theta grad = gradient(X, theta, Y) theta = theta - eta*grad if abs(cost_function(X, last_theta, Y) - cost_function(X, theta, Y)) < 1e-18: break return theta # theta = nomal(X, Y) # 闭合形式求解 theta = gradient_descent(X, theta, Y, 0.001) # 梯度下降 # 评估项(R2) def evaluation(x, theta, y, y_var): return 1 - ((np.sum((np.dot(x, theta)-y)**2))/(y_var*len(y))) MSE = np.sum(np.power((np.dot(x, theta)-y), 2))/len(y) cost = cost_function(x, theta, y) # 损失函数值 R2 = evaluation(x, theta, y, y_var) # 评估值 end = time.time() # print("The normal equations:") print("Gradient descent:") print("theta=") print(theta) print("MSE=", MSE) print("cost=", cost) print("R2=", R2) print('Running time: %s Seconds' % (end-start))
实验结果:
(一)正规方程求解
(二)梯度下降求解
可以看到两个方法得出的结果差别不大,用测试集进行测试时候,损失函数值均为0.18,评估项R²均为0.6多,梯度下降的拟合效果会比正规方程的好一点。在运算过程中,能很明显看到正规方程的计算速度要比梯度下降快很多,原因在于梯度下降在更新θ时候需要迭代很多次才能得到较优解。但是梯度下降在特征数量n较大时也能很好使用,而正规方程需要计算(X*X)-1,如果特征数量太多则运算代价较大因为矩阵的运算时间复杂度为O(n³),而且只适用于线性模型,不适用于逻辑回归模型等其他模型。在这个模型里面,由于特征数量不是很多,因此用正规式求解比较合理。