前言
上一篇我们对数据进行了读取并进行了可视化，今天我们来继续实现算法。
完整代码会在最后给出，如果你直接复制下面零散的代码可能会运行不了。
这篇的代码已经默认import了pandas,numpy等模块。
数据标准化
首先我们先取要用的三列数据作为训练数据集：
trainData = cReader[['mpg','displacement','acceleration']]
trainData.insert(0,'ones',1)
cols = trainData.shape[1]
X = trainData.iloc[:,0:cols-1]
Y = trainData.iloc[:,cols-1:cols]
X = np.mat(X.values)
Y = np.mat(Y.values)
for i in range(1,3):
X[:,i] = (X[:,i] - min(X[:,i])) / (max(X[:,i]) - min(X[:,I]))
Y[:,0] = (Y[:,0] - min(Y[:,0])) / (max(Y[:,0]) - min(Y[:,0]))
打印一下trainData前五行数据:
trainData数据
代码说明
在训练数据第一列加一列全为1的列矩阵trainData.insert(0,'ones',1)，目的是为了进行线性回归的时候简化计算，因为我们的目标方程是h(x)= θ0+ θ1𝓧1+ θ2𝓧2的一个常数项可以看成是θ0和1的乘积，其他项为θi和𝓧i的乘积
下图βi即为我们的θi
解释
cols得到trainData的列数，shape[0],shape[1]分别是trainData的行数和列数，下面把它的前三行'ones','mpg','displacement'分配给自变量矩阵𝓧最后一列分配给因变量矩阵Y
得到数据之后将它们标准化，关于做线性回归是否需要标准化数据，这里有一个比较好的解释
一般来说，我们再做线性回归时并不需要中心化和标准化数据。大多数情况下数据中的特征会以不同的测量单位展现，无论有没有中心化或者标准化都不会影响线性回归的结果。
因为估计出来的参数值β会恰当地将每个解释变量的单位x转化为响应变量的单位y.
但是标准化数据能将我们的结果更具有可解释性，比如β1=0.6
和β2=0.3, 我们可以理解为第一个特征的重要性是第二个特征的两倍。
这里采用以下公式进行标准化，对应上面代码的最后三行：
离差标准化公式
开始回归
首先用最小二乘法计算回归系数，并给出梯度下降的代价函数的代码表示
最小二乘法公式:
最小二乘法
代码:
theta_n = (X.T*X).I*X.T*Y
print(theta_n)
打印结果: [[ 0.58007057] [-0.0378746 ] [-0.35473364]]，从左到右分别是θ0，θ1，θ2
定义代价函数:
公式：
代价函数公式
代码:
我是直接自己定义了一个新模块linearRegrassion.py，在里面写了线性回归相关的函数，这也是上篇文章中开头的import linearRegrassion as lg
记得在新文件linearRegrassion.py里要
import pandas as pd
import numpy as np
代价函数costFunc：
def costFunc(X,Y,theta):
inner = np.power((X*theta.T)-Y,2)
return np.sum(inner)/(2*len(X))
观察公式，h(x)是预测值，y是实际值，通过他们俩的差来体现不同的拟合曲线的可靠程度，这个值越小说明对应的系数θ拟合出的曲线越接近实际，我们下面要做的就是用梯度下降的算法，取一系列θ值来逼近这个最优解，最后得到回归曲线。
梯度下降算法
公式:
梯度下降公式
解释
我们举一个代价函数的例子，看下面这个函数图像:
cost
我们在上面取一点(-30, 2500), 假设它对应一组我们待求的系数θ，此时代价函数值为2500，远没有达到最小的情况。
由于这个函数在(-∞, 20)区间内递减，所以我们肯定需要增大θ值来逼近最优解，于是我们对θ求偏导，然后再用θ减去偏导值，这样就可以使θ增大
(因为此时斜率k为负，所以减去之后肯定是增大的，相应的，如果这个点跑到对称轴右边，斜率变成正数，那么θ则会减少，仍然向最优解方向逼近)
说了这么多，我们取两个点画个图来看看:
k
梯度下降迭代到最后，斜率越来越接近0，代价函数也越来越接近最优解，此时我们得到的便是拟合度最高的系数θ
上两图的代码
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-100,+100,10)
y = pow(x-20,2)
k = 2*(x-20)
y1 = -20*(x-10)+100
y2 = -100*(x+30)+2500
plt.plot(x,y)
plt.plot(x,y1,label= 'k=-20')
plt.plot(x,y2,label= 'k=-100')
leg = plt.legend(loc='upper right',ncol=1)
plt.show()
梯度下降代码编写,还是在linearRegrassion.py模块中:
def gradientDescent(X,Y,theta,alpha,iters):
temp = np.mat(np.zeros(theta.shape))
cost = np.zeros(iters)
thetaNums = int(theta.shape[1])
print(thetaNums)
for i in range(iters):
error = (X*theta.T-Y)
for j in range(thetaNums):
derivativeInner = np.multiply(error,X[:,j])
temp[0,j] = theta[0,j] - (alpha*np.sum(derivativeInner)/len(X))
theta = temp
cost[i] = costFunc(X,Y,theta)
return theta,cost
解释一下几个参数:
alpha是公式中的𝒂，我们通常称之为步长，它决定了θ逼近最优解的速度, 过大会导致函数无法收敛得不到最优解，过小则会使逼近速度过慢。
iters是迭代次数，足够大的情况下得到的θ才有说服力
theta则是初始化迭代时给的一组θ值，最后得到拟合度最高的θ并在函数中返回
开始回归
开始回归之前，我们先设置一下学习的参数:
theta = np.mat([0,0,0])
iters = 100000
alpha = 0.001
回归结果及可视化
finalTheta,cost = lg.gradientDescent(X,Y,theta,alpha,iters)
print(finalTheta)
print(cost)
x1 = np.linspace(X[:,1].min(),X[:,1].max(),100)
x2 = np.linspace(X[:,2].min(),X[:,2].max(),100)
x1,x2 = np.meshgrid(x1,x2)
f = finalTheta[0,0] + finalTheta[0,1]*x1 + finalTheta[0,2]*x2
fig = plt.figure()
Ax = Axes3D(fig)
Ax.plot_surface(x1, x2, f, rstride=1, cstride=1, cmap=cm.viridis,label='prediction')
Ax.scatter(X[:100,1],X[:100,2],Y[:100,0],c='y')
Ax.scatter(X[100:250,1],X[100:250,2],Y[100:250,0],c='r')
Ax.scatter(X[250:,1],X[250:,2],Y[250:,0],c='b')
Ax.set_zlabel('acceleration') # 坐标轴
Ax.set_ylabel('displacement')
Ax.set_xlabel('mpg')
plt.show()
finalTheta，也就是最终的θ:
[[ 15.47017446 2.99065096 -3.31870705]]
最小二乘法估计的结果
[[ 17.74518558]
[ -0.63629322]
[ -5.95952521]]
代价函数的值:
[ 124.67279846 124.37076615 124.06949161 ..., 2.77449658 2.77449481
2.77449304]
可以看到迭代十万次之后得到的θ是和之前估计的比较接近的，如果你把迭代次数改为100万次，虽然计算时间长一点但可以看到结果是基本一样的
可视化:
回归结果
比较密集的部分还是基本分布在平面上下的
梯度下降的过程可视化:
代码：
fig, bx = plt.subplots(figsize=(8,6))
bx.plot(np.arange(iters), cost, 'r')
bx.set_xlabel('Iterations')
bx.set_ylabel('Cost')
bx.set_title('Error vs. Training Epoch')
plt.show()
梯度下降
如图，随着迭代次数的增加，代价函数越来越小，趋势越来越平稳，同时逼近最优解。
完整代码:
linearRegression.py
import pandas as pd
import numpy as np
def costFunc(X,Y,theta):
inner = np.power((X*theta.T)-Y,2)
return np.sum(inner)/(2*len(X))
def gradientDescent(X,Y,theta,alpha,iters):
temp = np.mat(np.zeros(theta.shape))
cost = np.zeros(iters)
thetaNums = int(theta.shape[1])
print(thetaNums)
for i in range(iters):
error = (X*theta.T-Y)
for j in range(thetaNums):
derivativeInner = np.multiply(error,X[:,j])
temp[0,j] = theta[0,j] - (alpha*np.sum(derivativeInner)/len(X))
theta = temp
cost[i] = costFunc(X,Y,theta)
return theta,cost
lgScript.py
from io import StringIO
from urllib import request
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
import ssl
import pandas as pd
import numpy as np
import linearRegrassion as lg
ssl._create_default_https_context = ssl._create_unverified_context
names =["mpg","cylinders","displacement","horsepower",
"weight","acceleration","model year","origin","car name"]
url = 'https://archive.ics.uci.edu/ml/machine-learning-databases/auto-mpg/auto-mpg.data'
s = request.urlopen(url).read().decode('utf8')
dataFile = StringIO(s)
cReader = pd.read_csv(dataFile,delim_whitespace=True,names=names)
ax = plt.subplot(111, projection='3d') # 创建一个三维的绘图工程
ax.scatter(cReader["mpg"][:100],cReader["displacement"][:100],cReader["acceleration"][:100],c='y')
ax.scatter(cReader["mpg"][100:250],cReader["displacement"][100:250],cReader["acceleration"][100:250],c='r')
ax.scatter(cReader["mpg"][250:],cReader["displacement"][250:],cReader["acceleration"][250:],c='b')
ax.set_zlabel('acceleration') # 坐标轴
ax.set_ylabel('displacement')
ax.set_xlabel('mpg')
plt.show()
plt.scatter(cReader["mpg"],cReader["displacement"])
plt.xlabel('mpg')
plt.ylabel('displacement')
plt.show()
trainData = cReader[['mpg','displacement','acceleration']]
trainData.insert(0,'ones',1)
print(trainData.head(5))
cols = trainData.shape[1]
X = trainData.iloc[:,0:cols-1]
Y = trainData.iloc[:,cols-1:cols]
X = np.mat(X.values)
Y = np.mat(Y.values)
for i in range(1,3):
X[:,i] = (X[:,i] - min(X[:,i])) / (max(X[:,i]) - min(X[:,i]))
print(X[:5:,:3])
#Y[:,0] = (Y[:,0] - min(Y[:,0])) / (max(Y[:,0]) - min(Y[:,0]))
print(Y[:5,0])
theta_n = (X.T*X).I*X.T*Y
print(theta_n)
theta = np.mat([0,0,0])
iters = 100000
alpha = 0.001
finalTheta,cost = lg.gradientDescent(X,Y,theta,alpha,iters)
print(finalTheta)
print(cost)
x1 = np.linspace(X[:,1].min(),X[:,1].max(),100)
x2 = np.linspace(X[:,2].min(),X[:,2].max(),100)
x1,x2 = np.meshgrid(x1,x2)
f = finalTheta[0,0] + finalTheta[0,1]*x1 + finalTheta[0,2]*x2
fig = plt.figure()
Ax = Axes3D(fig)
Ax.plot_surface(x1, x2, f, rstride=1, cstride=1, cmap=cm.viridis,label='prediction')
Ax.scatter(X[:100,1],X[:100,2],Y[:100,0],c='y')
Ax.scatter(X[100:250,1],X[100:250,2],Y[100:250,0],c='r')
Ax.scatter(X[250:,1],X[250:,2],Y[250:,0],c='b')
Ax.set_zlabel('acceleration') # 坐标轴
Ax.set_ylabel('displacement')
Ax.set_xlabel('mpg')
plt.show()
fig, bx = plt.subplots(figsize=(8,6))
bx.plot(np.arange(iters), cost, 'r')
bx.set_xlabel('Iterations')
bx.set_ylabel('Cost')
bx.set_title('Error vs. Training Epoch')
plt.show()
总结
这个系列算是学习吴恩达机器学习的笔记与作业，由于平时课程较多，所以不定期更新，下一篇大概是用Python实现逻辑回归，希望不足之处大家可以讨论一下。
如果觉得《python求解多元方程最优解_Python实现梯度下降算法求多元线性回归(二)》对你有帮助，请点赞、收藏，并留下你的观点哦！
