线性回归代码实现

引入

上一节讲解了线性回归的基本原理，这一节主要通过代码加深学习和理解，本节代码是在python3.7集成的IDLE上完成~

代码和数据放在百度云盘（永久有效）

 链接:https://pan.baidu.com/s/1Mb4nKGRvHHIXGe95wzx_Ag 

提取码:a1zt

---

加载数据

# 导包
import numpy as np
import matplotlib.pyplot as plt

# 定义加载数据函数
def loadData():
    data = np.loadtxt('dataLinearRegression.txt',delimiter=',')
    n = data.shape[1]-1 # 特征数
    X = data[:,0:n]
    y = data[:,-1].reshape(-1,1)
    return X,y

# 测试
# >>> import LinearRegression as lr
# >>> X,y = lr.loadData()
# >>> X.shape # X是一维
# (97, 1) 
# >>> y.shape # y是一维
# (97, 1)
# >>> import matplotlib.pyplot as plt
# >>> plt.scatter(X,y)
# <matplotlib.collections.PathCollection object at 0x000001BDE10D9E08>
# >>> plt.show()

数据如下图所示：

---

最小二乘法代码

# 最小二乘法
X_origin,y = loadData()
# 还记得上一节为了计算方便，我们加入x0=1
X = np.insert(X_origin,0,values=1,axis=1)

# 在上一节中，使用最小二乘法得到theta的计算公式，忘记的记得去复习

theta = np.linalg.inv(X.T.dot(X)).dot(X.T).dot(y)
# numpy.linalg模块包含线性代数的函数，使用这个模块，可以计算逆矩阵、特征值、行列式等
# np.linalg.inv()# 求逆矩阵

# 自行测试    
# >>> import LinearRegression as lr
# >>> lr.theta
# array([[-3.89578088],
#        [ 1.19303364]])

# 至此，参数theta求出，即得到映射关系

import matplotlib.pyplot as plt
plt.scatter(X_origin,y)
h_theta = theta[0]+theta[1]*X_origin
plt.plot(X_origin,h_theta)
plt.show()

# 自行测试
#>>> import LinearRegression as lr

可以看出拟合的效果还是很不错的~

---

梯度下降法代码

# 梯度下降法
# 梯度下降法是对损失函数不断迭代，求出损失函数极值的问题

# 定义梯度下降法
def gradientDescent(X,y,theta,iterations,alpha):
    c = np.ones(X.shape[0]).transpose() # 插入x0=1
    X = np.insert(X,0,values=c,axis=1) # 对原始数据加入一个全为1的列
    m = X.shape[0] # 数据条数
    n = X.shape[1] # 特征数
    costs = np.zeros(iterations)
    for num in range(iterations):
        for j in range(n):
            theta[j] = theta[j]+(alpha/m)*(np.sum((y-np.dot(X,theta))*X[:,j].reshape(-1,1))) # 比对公式，认真理解
    return theta

X,y = loadData()
theta = np.zeros(X.shape[1]+1).reshape(-1,1) # 这里加1，是因为加入x0=1这一列，相当于特征多了一个
iterations = 400
alpha = 0.01
theta = gradientDescent(X,y,theta,iterations,alpha)

# 画拟合直线和数据散点图
plt.scatter(X,y)
h_theta = theta[0]+theta[1]*X # 最终的预测函数
plt.plot(X,h_theta)
plt.show()

# 自行测试
# >>> import LinearRegression as lr

拟合的效果还是可以的~

下面我们看看损失函数的变化趋势

# 定义损失函数
def computeCost(X,y,theta):
    m = X.shape[0] # 有多少行数据
    return np.sum(np.power(np.dot(X,theta)-y,2))/(2*m)

# 定义梯度下降法
def gradientDescent(X,y,theta,iterations,alpha):
    c = np.ones(X.shape[0]).transpose() # 插入x0=1
    X = np.insert(X,0,values=c,axis=1) # 对原始数据加入一个全为1的列
    m = X.shape[0] # 数据条数
    n = X.shape[1] # 特征数
    costs = np.zeros(iterations)
    for num in range(iterations):
        for j in range(n):
            theta[j] = theta[j]+(alpha/m)*(np.sum((y-np.dot(X,theta))*X[:,j].reshape(-1,1)))
        costs[num] = computeCost(X,y,theta)
    return theta,costs
X,y = loadData()
theta = np.zeros(X.shape[1]+1).reshape(-1,1)
iterations = 400
alpha = 0.01
theta,costs = gradientDescent(X,y,theta,iterations,alpha)

# 画损失函数图
x_axis = np.linspace(1,iterations,iterations) # 用来创建等差数列，x轴从1到400，样本数量400个
plt.plot(x_axis,costs[0:iterations])
plt.show()

# 自行测试
# >>> import LinearRegression as lr

可能数据分布的原因，导致损失函数下降过程没有那么平滑

下面对数据进行一下归一化

# 归一化
def featureNormalize(X):
    mu = np.average(X,axis=0)# 均值
    sigma = np.std(X,axis=0,ddof=1)# 方差
    X = (X-mu)/sigma
    return X,mu,sigma

# 定义损失函数
def computeCost(X,y,theta):
    m = X.shape[0] # 有多少行数据
    return np.sum(np.power(np.dot(X,theta)-y,2))/(2*m)

# 定义梯度下降法
def gradientDescent(X,y,theta,iterations,alpha):
    c = np.ones(X.shape[0]).transpose() # 插入x0=1
    X = np.insert(X,0,values=c,axis=1) # 对原始数据加入一个全为1的列
    m = X.shape[0] # 数据条数
    n = X.shape[1] # 特征数
    costs = np.zeros(iterations)
    for num in range(iterations):
        for j in range(n):
            theta[j] = theta[j]+(alpha/m)*(np.sum((y-np.dot(X,theta))*X[:,j].reshape(-1,1)))
        costs[num] = computeCost(X,y,theta)
    return theta,costs

X_origin,y = loadData() # 加载数据
X,mu,sigma = featureNormalize(X_origin) # 归一化
theta = np.zeros(X.shape[1]+1).reshape(-1,1) # 初始化theta
iterations = 400
alpha = 0.01
theta,costs = gradientDescent(X,y,theta,iterations,alpha)

# 画拟合直线和数据散点图
plt.scatter(X,y)
h_theta = theta[0]+theta[1]*X
plt.plot(X,h_theta)
plt.show()

# 自行测试
# >>> import LinearRegression as lr

# 画损失函数图
x_axis = np.linspace(1,iterations,iterations) # 用来创建等差数列，x轴从1到400，样本数量400个
plt.plot(x_axis,costs[0:iterations])
plt.show()

# 自行测试
# >>> import LinearRegression as lr

可以看出损失函数平稳减少~

如果一个新的数据到来，我们如何预测呢？

# 预测函数
def predict(X):
    X = (X-mu)/sigma
    c = np.ones(X.shape[0]).transpose()
    X = np.insert(X,0,values=c,axis=1)
    return np.dot(X,theta)

# >>> import LinearRegression as lr
# >>> print(lr.predict([[6.1101]]))
# [[3.3347208]]

---

总结

本节主要是对最小二乘法和梯度下降法进行代码实现，需要注意的是主要当数据分布差异较大时，需要对数据进行归一化。

欢迎交流~