一:梯度下降法推导公式:
二:当训练集为1维时
#进行数据分析所需库,可以看做是对numpy工具的补充
import pandas as pd
import numpy as np
#应该把Seaborn视为matplotlib的补充,作图所用工具,在大多数情况下使用seaborn就能做出很具有吸引力的图,而使用matplotlib就能制作具有更多特色的图
import seaborn as sns
import matplotlib.pyplot as plt
#设置绘画的图标格式和颜色
sns.set(context="notebook", style="whitegrid", palette="dark")
#读取数据并赋予列名,此时df共有两列,列名分别为population和profit
df = pd.read_csv('ex1data1.txt', names=['population', 'profit'])
def get_x(df):
#ones:格式化为一个m列的全为1的向量
ones = pd.DataFrame({'ones':np.ones(len(df))})
#将date格式化为ones向量右边加上df矩阵的矩阵
data = pd.concat([ones, df], axis=1)
#返回data矩阵的前两列,返回m * 2的矩阵
return data.iloc[:, :-1].as_matrix() # 这个操作返回 ndarray,不是矩阵
def get_y(df):
#返回df矩阵的最后一列,返回m × 1的向量
return np.array(df.iloc[:, -1])
def normalize_feature(df):
# 特征缩放,对df中的两列数据分别进行特征缩放 :(x - x平均值)/ x方差
return df.apply(lambda column: (column - column.mean())/column.std())
X = get_x(df)
print(X.shape, type(X)) #X为m × 2的矩阵
y = get_y(df)
print(y.shape, type(y)) #y为m × 1的向量
theta = np.zeros(X.shape[1])#X.shape[1]=2,代表特征数n X.shape(0) = 97;X.shape(1) = 2;将theta初始化为2行的零向量(列向量)
def lr_cost(theta, X, y):
# """
# 计算theta固定时,此时的代价函数值
# X: R(m*n), m 样本数, n 特征数
# y: R(m)
# theta : R(n), 线性回归的参数
# """
m = X.shape[0]#m为样本数
#X矩阵和theta矩阵的点积所得矩阵(m × 1) - 矩阵y(m × 1)
inner = X @ theta - y # R(m*1),X @ theta等价于X.dot(theta)
#square_sum为inner的每个元素平方之和,即二范数
square_sum = inner.T @ inner
#得到此时代价函数的值
cost = square_sum / (2 * m)
return cost
#将代价函数对theta求导,即求梯度
def gradient(theta, X, y):
#m为样本数
m = X.shape[0]
#这里的inner为代价函数对theta求导的结果
inner = X.T @ (X @ theta - y) # (m,n).T @ (m, 1) -> (n, 1),X @ theta等价于X.dot(theta)
return inner / m
#梯度下降函数
def batch_gradient_decent(theta, X, y, epoch, alpha=0.01):
# 拟合线性回归,返回参数和代价
# epoch: 批处理的轮数
# alpha: theta移动的步长
# """
#得到初始theta时的代价
cost_data = [lr_cost(theta, X, y)]
# 拷贝一份,不和原来的theta混淆
_theta = theta.copy()
#开始迭代epoch次
for _ in range(epoch):
#根据当前梯度更新theta
_theta = _theta - alpha * gradient(_theta, X, y)
#记录此时的代价
cost_data.append(lr_cost(_theta, X, y))
return _theta, cost_data
#迭代500次求最小代价和所对应的theta
epoch = 500
final_theta, cost_data = batch_gradient_decent(theta, X, y, epoch)
# 计算最终的代价
lr_cost(final_theta, X, y)
#画出代价函数值变化图
#可以看到从第二轮代价数据变换很大,接下来平稳了
ax = sns.tsplot(cost_data, time=np.arange(epoch+1))
ax.set_xlabel('epoch')
ax.set_ylabel('cost')
plt.show()
b = final_theta[0] # intercept,Y轴上的截距
m = final_theta[1] # slope,斜率
#画出原数据点和线性回归的最终结果
plt.scatter(df.population, df.profit, label="Training data")
plt.plot(df.population, df.population*m + b, label="Prediction")
plt.legend(loc=2)
plt.show()
三:当训练集为多维时
当训练集为多维时只需将上述代码稍作改动,完整代码如下:
#进行数据分析所需库,可以看做是对numpy工具的补充
import pandas as pd
import numpy as np
#应该把Seaborn视为matplotlib的补充,作图所用工具,在大多数情况下使用seaborn就能做出很具有吸引力的图,而使用matplotlib就能制作具有更多特色的图
import seaborn as sns
import matplotlib.pyplot as plt
#设置绘画的图标格式和颜色
sns.set(context="notebook", style="whitegrid", palette="dark")
#读取数据并赋予列名,此时df共有两列,列名分别为population和profit
df = pd.read_csv('ex1data2.txt', names=['square', 'bedrooms', 'price'])
df.head()
def get_x(df):
#ones:格式化为一个m列的全为1的向量
ones = pd.DataFrame({'ones':np.ones(len(df))})
#将date格式化为ones向量右边加上df矩阵的矩阵
data = pd.concat([ones, df], axis=1)
#返回data矩阵的前两列,返回m * 2的矩阵
return data.iloc[:, :-1].as_matrix() # 这个操作返回 ndarray,不是矩阵
def get_y(df):
#返回df矩阵的最后一列,返回m × 1的向量
return np.array(df.iloc[:, -1])
def normalize_feature(df):
# 特征缩放,对df中的两列数据分别进行特征缩放 :(x - x平均值)/ x方差
return df.apply(lambda column: (column - column.mean())/column.std())
data = normalize_feature(df) #特征缩放
X = get_x(data)
print(X.shape, type(X)) #X为m × 2的矩阵
y = get_y(data)
print(y.shape, type(y)) #y为m × 1的向量
theta = np.zeros(X.shape[1])#X.shape[1]=2,代表特征数n X.shape(0) = 97;X.shape(1) = 2;将theta初始化为2行的零向量(列向量)
def lr_cost(theta, X, y):
#
# 计算theta固定时,此时的代价函数值
# X: R(m*n), m 样本数, n 特征数
# y: R(m)
# theta : R(n), 线性回归的参数
#
m = X.shape[0]#m为样本数
#X矩阵和theta矩阵的点积所得矩阵(m × 1) - 矩阵y(m × 1)
inner = X @ theta - y # R(m*1),X @ theta等价于X.dot(theta)
#square_sum为inner的每个元素平方之和,即二范数
square_sum = inner.T @ inner
#得到此时代价函数的值
cost = square_sum / (2 * m)
return cost
#将代价函数对theta求导,即求梯度
def gradient(theta, X, y):
#m为样本数
m = X.shape[0]
#这里的inner为代价函数对theta求导的结果
inner = X.T @ (X @ theta - y) # (m,n).T @ (m, 1) -> (n, 1),X @ theta等价于X.dot(theta)
return inner / m
#梯度下降函数
def batch_gradient_decent(theta, X, y, epoch, alpha=0.01):
# 拟合线性回归,返回参数和代价
# epoch: 批处理的轮数
# alpha: theta移动的步长
#
#得到初始theta时的代价
cost_data = [lr_cost(theta, X, y)]
# 拷贝一份,不和原来的theta混淆
_theta = theta.copy()
#开始迭代epoch次
for _ in range(epoch):
#根据当前梯度更新theta
_theta = _theta - alpha * gradient(_theta, X, y)
#记录此时的代价
cost_data.append(lr_cost(_theta, X, y))
return _theta, cost_data
base = np.logspace(-1, -5, num=4)
candidate = np.sort(np.concatenate((base, base*3)))
epoch = 50
fig, ax = plt.subplots(figsize=(16, 9))
for alpha in candidate:
_, cost_data = batch_gradient_decent(theta, X, y, epoch, alpha=alpha)
ax.plot(np.arange(epoch+1), cost_data, label=alpha)
ax.set_xlabel('epoch', fontsize=18)
ax.set_ylabel('cost', fontsize=18)
ax.legend(bbox_to_anchor=(1.05, 1), loc=2, borderaxespad=0.)
ax.set_title('learning rate', fontsize=18)
plt.show()
'''
#画出代价函数值变化图
#可以看到从第二轮代价数据变换很大,接下来平稳了
ax = sns.tsplot(cost_data, time=np.arange(epoch+1))
ax.set_xlabel('epoch')
ax.set_ylabel('cost')
plt.show()
'''