参考：某大佬的github
建议打开上面的链接看看。
此篇主要补充一些公式转化成线性代数语言和实现。
如有纰漏，欢迎指正。
主要部分有：加载数据，归一化数据，计算代价cost，梯度下降。
另外需要一个训练函数，一个预测函数，一些可视化函数和主函数。
x-特征
y-真实结果
h-假设函数的结果（计算得到的结果）
m-训练样本的数量
n-特征的维数（比如 2 ：房屋面积，使用年限）
theta-参数
去上面的链接复制粘贴就行
def load_data(filename, delimiter, dtype):
loaddata = np.loadtxt(filename, delimiter=delimiter, dtype=dtype)
# print(loaddata)
return loaddata
归一化数据
公式：
上面的链接使用的标准差（代码std)，而我用的最大值-最小值（代码ptp)
def normalization(data_x): # 注意：接收的是x 此函数把参数转化成-0.5到+0.5之间
# 使用x = （x-平均值）/范围 替代x ---------- 也可以除以标准差
data_x_norm = np.array(data_x) # 将X转化为numpy数组对象，才可以进行矩阵的运算
mean = np.mean(data_x_norm, 0) # 求每一列的平均值（0指定为列，1代表行）（一列是同一种特征）
# print(f'mean={mean}')
ptp = np.ptp(data_x_norm, 0) # 求每一列的范围
# print(f'ptp={ptp}')
for col in range(data_x_norm.shape[1]): # python遍历np数组的列
# shape[0]第一维度的长度相当于行数 shape[1]第二维度的长度相当于列数
data_x_norm[:, col] = (data_x_norm[:, col]-mean[col])/ptp[col]
# 其中data_x_norm[:, col]应该是一列，mean[col]，ptp[col] 应该是个数字
# print(f'data_x_norm[:, col] ={data_x_norm[:, col] }') # 显示的应该是科学计数法表示的数据，-0.5到+0.5之间是正常的
return data_x_norm, mean, ptp
计算代价cost
有几点要注意的：
1.特征x应该有一列1，实现如下：
x_fin = np.hstack((np.ones((m, 1)), x_norm)) # 在X前加一列1 原因见吴恩达课 # np.vstack():在竖直方向上堆叠 np.hstack():在水平方向上平铺
原因吴恩达老师讲过，这样做可以使得 输出h = x 乘 参数theta
计算输出（假设函数）h的方法见上图，代码如下：
h = np.dot(x, theta) # 计算预测值（把数据代入假设）
计算代价cost的公式：
对应代码：
cost_j = np.transpose(np.dot(x_cost, theta_cost) - y_cost) * (np.dot(x_cost, theta_cost) - y_cost) / (2 * m)
其中h我们刚刚讲过计算方法，
平方和的计算需要应用线性代数知识：a转置乘a就是平方和
证明如下（2022李永乐线代33页）
计算代价cost的代码如下：
def cost(x_cost, y_cost, m, theta_cost): # 注意x第一列应该全为1
# print(f'x_cost={x_cost}')
# print(f'y_cost={y_cost}')
# print(f'theta_cost={theta_cost}')
cost_j = np.transpose(np.dot(x_cost, theta_cost) - y_cost) * (np.dot(x_cost, theta_cost) - y_cost) / (2 * m)
# 计算代价J # at*a等价于平方和 线代知识
return cost_j # 应该是个行向量？
梯度下降
求导可以自己试试，很简单的复合函数求导
公式中
这一部分的处理其实和刚刚很像，h还用刚刚的方法计算：
h = np.dot(x, theta) # 计算预测值（把数据代入假设）
剩下的“乘积求和” 和 “平方和”是一样的，用（h-y）转置乘x就是“乘积求和”的结果了。
实现如下：
def gradient_descent(x, y, m, n, theta, alpha, num_iters):
temp = np.matrix(np.zeros((n, num_iters))) # 暂存每次迭代计算的theta，转化为矩阵形式
j_history = np.zeros((num_iters, 1)) # 记录每次迭代计算的代价值
for i in range(num_iters): # 遍历迭代次数
h = np.dot(x, theta) # 计算预测值（把数据代入假设）
temp[:, i] = theta - ((alpha / m) * (np.dot(np.transpose(x), h - y))) # 梯度的计算
# temp[:, i]遍历列 transpose(x) 取转置就可以一一对应相乘了
theta = temp[:, i]
j_history[i] = cost(x, y, m, theta) # 调用计算代价函数
print('.', end=' ')
print()
return theta, j_history
完整代码
import numpy as np
from matplotlib import pyplot as plt
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:windowsfontssimsun.ttc", size=14) # 解决windows环境下画图汉字乱码问题
def load_data(filename, delimiter, dtype):
loaddata = np.loadtxt(filename, delimiter=delimiter, dtype=dtype)
# print(loaddata)
return loaddata
def normalization(data_x): # 注意：接收的是x 此函数把参数转化成-0.5到+0.5之间
# 使用x = （x-平均值）/范围 替代x ---------- 也可以除以标准差
data_x_norm = np.array(data_x) # 将X转化为numpy数组对象，才可以进行矩阵的运算
mean = np.mean(data_x_norm, 0) # 求每一列的平均值（0指定为列，1代表行）（一列是同一种特征）
# print(f'mean={mean}')
ptp = np.ptp(data_x_norm, 0) # 求每一列的范围
# print(f'ptp={ptp}')
for col in range(data_x_norm.shape[1]): # python遍历np数组的列
# shape[0]第一维度的长度相当于行数 shape[1]第二维度的长度相当于列数
data_x_norm[:, col] = (data_x_norm[:, col]-mean[col])/ptp[col]
# 其中data_x_norm[:, col]应该是一列，mean[col]，ptp[col] 应该是个数字
# print(f'data_x_norm[:, col] ={data_x_norm[:, col] }') # 显示的应该是科学计数法表示的数据，-0.5到+0.5之间是正常的
return data_x_norm, mean, ptp
def cost(x_cost, y_cost, m, theta_cost): # 注意x第一列应该全为1
# print(f'x_cost={x_cost}')
# print(f'y_cost={y_cost}')
# print(f'theta_cost={theta_cost}')
cost_j = np.transpose(np.dot(x_cost, theta_cost) - y_cost) * (np.dot(x_cost, theta_cost) - y_cost) / (2 * m)
# 计算代价J # at*a等价于平方和 线代知识
return cost_j # 应该是个行向量？
def gradient_descent(x, y, m, n, theta, alpha, num_iters):
temp = np.matrix(np.zeros((n, num_iters))) # 暂存每次迭代计算的theta，转化为矩阵形式
j_history = np.zeros((num_iters, 1)) # 记录每次迭代计算的代价值
for i in range(num_iters): # 遍历迭代次数
h = np.dot(x, theta) # 计算预测值（把数据代入假设）
temp[:, i] = theta - ((alpha / m) * (np.dot(np.transpose(x), h - y))) # 梯度的计算
# temp[:, i]遍历列 transpose(x) 取转置就可以一一对应相乘了
theta = temp[:, i]
j_history[i] = cost(x, y, m, theta) # 调用计算代价函数
print('.', end=' ')
print()
return theta, j_history
def linear_regression(alpha=0.01, num_iters=400): # 训练
data = load_data('data.csv', ',', np.float64) # 加载数据
# 拆解数据
x = data[:, 0:-1] # 3列，0，1是【x】，2是【y】
y = data[:, -1] # y对应最后一列
x_norm, mean, ptp = normalization(x) # 归一化x
plot_x1_x2(x_norm) # 画图
y = y.reshape(-1, 1) # 将行向量转化为列
m = len(y) # 训练样本的数量
x_fin = np.hstack((np.ones((m, 1)), x_norm)) # 在X前加一列1 原因见吴恩达课 # np.vstack():在竖直方向上堆叠 np.hstack():在水平方向上平铺
n = data.shape[1] # data的列数 即特征的维数（2+1 有一列1）
theta = np.zeros((n, 1)) # 初始化参数
print(u"n执行梯度下降算法....n")
theta, j_history = gradient_descent(x_fin, y, m, n, theta, alpha, num_iters)
plotJ(j_history, num_iters)
return theta, mean, ptp # 返回学习的结果theta 平均值mean 范围ptp(用于预测时用normalization(x)同样的放缩处理输入数据)
def plotJ(j_history, num_iters): # 画每次迭代代价的变化图
x = np.arange(1, num_iters+1)
plt.plot(x, j_history)
plt.xlabel(u"迭代次数", fontproperties=font) # 注意指定字体，要不然出现乱码问题
plt.ylabel(u"代价值", fontproperties=font)
plt.title(u"代价随迭代次数的变化", fontproperties=font)
plt.show()
def plot_x1_x2(x):
plt.scatter(x[:, 0], x[:, 1])
plt.show()
def predict(pre_x, theta, mean, ptp):
norm_pre_x = (pre_x - mean) / ptp
fin_pre_x = np.hstack((np.ones((1)), norm_pre_x))
result = np.dot(fin_pre_x, theta)
return result
if __name__ == "__main__":
theta, mean, ptp = linear_regression()
pre_x = np.array([1650, 4])
print(predict(pre_x, theta, mean, ptp))
