利用Python编程，分别使用梯度下降法和最小二乘法求解多元函数
使用Python编程，分别根据梯度法和最小二乘法求解多元函数问题
- 分别使用梯度下降法和最小二乘法求解多元函数并进行比较，这里使用jupyter notebook平台进行Python编程
- 一、题目描述
- 二、使用梯度下降法求解多元函数
- （一）梯度下降法基本原理
- （二）梯度下降法公式
- （三）Python代码实现
- 1、导入要使用到的库、定义变量并赋值
- 2、代价函数
- 3、使用梯度下降法求解多元函数系数
- 4、打印系数和方程
- 5、绘制线性回归方程拟合图
- 6、完整代码
- 7、运行结果
- 三、使用最小二乘法求解多元函数
- （一）最小二乘法基本原理
- （二）Python代码实现
- 1、自行推导过程编写代码
- 运行结果
- 2、采用sklearn库
- 运行结果
- （四）将梯度下降法与最小二乘法进行对比
分别使用梯度下降法和最小二乘法求解多元函数并进行比较，这里使用jupyter notebook平台进行Python编程
一、题目描述
为求得某个地区的商品店的月营业额与店铺的面积、该店距离车站距离的相关性大小，以店铺面积、距离车站的距离、以及月营业额建立线性回归方程，并求解方程，得到相关系数。
将表中的数据录入到Excel中，编程时会用到这些数据，如下：
二、使用梯度下降法求解多元函数
（一）梯度下降法基本原理
梯度下降法又称最速下降法，是求解无约束最优化问题的一种最常用的方法，在对损失函数最小化时经常使用。梯度下降法是一种迭代算法。选取适当的初值x(0)，不断迭代，更新x的值，进行目标函数的极小化，直到收敛。由于负梯度方向时使函数值下降最快的方向，在迭代的每一步，以负梯度方向更新x的值，从而达到减少函数值的目的。
（二）梯度下降法公式
（三）Python代码实现
1、导入要使用到的库、定义变量并赋值
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=np.genfromtxt('D:/area_distance.csv',delimiter=',')
x_data=data[:,:-1]
y_data=data[:,2]
#定义学习率、斜率、截据
#设方程为y=theta1x1+theta2x2+theta0
lr=0.00001
theta0=0
theta1=0
theta2=0
#定义最大迭代次数，因为梯度下降法是在不断迭代更新k与b
epochs=10000
2、代价函数
def compute_error(theta0,theta1,theta2,x_data,y_data):totalerror=0for i in range(0,len(x_data)):#定义一共有多少样本点totalerror=totalerror+(y_data[i]-(theta1*x_data[i,0]+theta2*x_data[i,1]+theta0))**2return totalerror/float(len(x_data))/2
3、使用梯度下降法求解多元函数系数
def gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs):m=len(x_data)for i in range(epochs):theta0_grad=0theta1_grad=0theta2_grad=0for j in range(0,m):theta0_grad-=(1/m)*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta2)+y_data[j])theta1_grad-=(1/m)*x_data[j,0]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])theta2_grad-=(1/m)*x_data[j,1]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])theta0=theta0-lr*theta0_gradtheta1=theta1-lr*theta1_gradtheta2=theta2-lr*theta2_gradreturn theta0,theta1,theta2
4、打印系数和方程
theta0,theta1,theta2=gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs)
print('迭代次数：{0} 学习率：{1}之后 a0={2},a1={3},a2={4},代价函数为{5}'.format(epochs,lr,theta0,theta1,theta2,compute_error(theta0,theta1,theta2,x_data,y_data)))
print("多元线性回归方程为:y=",theta1,"X1+",theta2,"X2+",theta0)
5、绘制线性回归方程拟合图
#画图
ax=plt.figure().add_subplot(111,projection='3d')
ax.scatter(x_data[:,0],x_data[:,1],y_data,c='r',marker='o')
x0=x_data[:,0]
x1=x_data[:,1]
#生成网格矩阵
x0,x1=np.meshgrid(x0,x1)
z=theta0+theta1*x0+theta2*x1
#画3d图
ax.plot_surface(x0,x1,z)
ax.set_xlabel('area')
ax.set_ylabel('distance')
ax.set_zlabel("Monthly turnover")
plt.show()
6、完整代码
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=np.genfromtxt('D:/area_distance.csv',delimiter=',')
x_data=data[:,:-1]
y_data=data[:,2]
#定义学习率、斜率、截据
#设方程为y=theta1x1+theta2x2+theta0
lr=0.00001
theta0=0
theta1=0
theta2=0
#定义最大迭代次数，因为梯度下降法是在不断迭代更新k与b
epochs=10000
#定义最小二乘法函数-损失函数（代价函数）
def compute_error(theta0,theta1,theta2,x_data,y_data):totalerror=0for i in range(0,len(x_data)):#定义一共有多少样本点totalerror=totalerror+(y_data[i]-(theta1*x_data[i,0]+theta2*x_data[i,1]+theta0))**2return totalerror/float(len(x_data))/2
#梯度下降算法求解参数
def gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs):m=len(x_data)for i in range(epochs):theta0_grad=0theta1_grad=0theta2_grad=0for j in range(0,m):theta0_grad-=(1/m)*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta2)+y_data[j])theta1_grad-=(1/m)*x_data[j,0]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])theta2_grad-=(1/m)*x_data[j,1]*(-(theta1*x_data[j,0]+theta2*x_data[j,1]+theta0)+y_data[j])theta0=theta0-lr*theta0_gradtheta1=theta1-lr*theta1_gradtheta2=theta2-lr*theta2_gradreturn theta0,theta1,theta2
#进行迭代求解
theta0,theta1,theta2=gradient_descent_runner(x_data,y_data,theta0,theta1,theta2,lr,epochs)
print('迭代次数：{0} 学习率：{1}之后 a0={2},a1={3},a2={4},代价函数为{5}'.format(epochs,lr,theta0,theta1,theta2,compute_error(theta0,theta1,theta2,x_data,y_data)))
print("多元线性回归方程为:y=",theta1,"X1+",theta2,"X2+",theta0)
#画图
ax=plt.figure().add_subplot(111,projection='3d')
ax.scatter(x_data[:,0],x_data[:,1],y_data,c='r',marker='o')
x0=x_data[:,0]
x1=x_data[:,1]
#生成网格矩阵
x0,x1=np.meshgrid(x0,x1)
z=theta0+theta1*x0+theta2*x1
#画3d图
ax.plot_surface(x0,x1,z)
ax.set_xlabel('area')
ax.set_ylabel('distance')
ax.set_zlabel("Monthly turnover")
plt.show()
7、运行结果
三、使用最小二乘法求解多元函数
（一）最小二乘法基本原理
在我们研究两个变量（x,y）之间的相互关系时，通常可以得到一系列成对的数据（x1,y1.x2,y2… xm,ym）；将这些数据描绘在x -y直角坐标系中，若发现这些点在一条直线附近，可以令这条直线方程为（式1-1）
其中：a0、a1 是任意实数
为建立这直线方程就要确定a0和a1，应用《最小二乘法原理》，将实测值Yi与利用计算值Yj（Yj=a0+a1Xi）（式1-1）的离差（Yi-Yj）的平方和 最小为“优化判据”。
令：φ = （式1-2)
把（式1-1）代入（式1-2）中得：
φ = （式1-3)
当 最小时，可用函数 φ 对a0、a1求偏导数，令这两个偏导数等于零。
∑2(a0 + a1Xi - Yi）=0（式1-4)
∑2Xi（a0 +a1Xi - Yi）=0（式1-5)
亦即：
na0 + （∑Xi ) a1 = ∑Yi （式1-6)
（∑Xi ) a0 + （∑Xi^2 ) a1 = ∑（Xi*Yi) （式1-7)
得到的两个关于a0、 a1为未知数的两个方程组，解这两个方程组得出：
a0 = （∑Yi) / n - a1（∑Xi) / n （式1-8)
a1 = [n∑(Xi Yi) - （∑Xi ∑Yi)] / (n∑Xi^2 -∑Xi∑Xi)（式1-9)
这时把a0、a1代入（式1-1）中， 此时的(式1-1）就是我们回归的一元线性方程即：数学模型。
在回归过程中，回归的关联式不可能全部通过每个回归数据点（x1,y1. x2,y2…xm,ym），为了判断关联式的好坏，可借助相关系数“R”，统计量“F”，剩余标准偏差“S”进行判断；“R”越趋近于 1 越好；“F”的绝对值越大越好；“S”越趋近于 0 越好。
R = [∑XiYi - m （∑Xi / m）（∑Yi / m)]/ SQR{[∑Xi2 - m （∑Xi / m)2][∑Yi2 - m （∑Yi / m)2]} （式1-10) *
在（式1-10）中，m为样本容量，即实验次数；Xi、Yi分别为任意一组实验数据X、Y的数值。
（二）Python代码实现
1、自行推导过程编写代码
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import seaborn as sns%matplotlib inline
data = np.genfromtxt("D:/area_distance.csv",delimiter=",")
X1=data[0:10,0]#面积
X2=data[0:10,1]#离车站的距离
Y=data[0:10,2]#月营业额
#将因变量赋值给矩阵Y1
Y1=np.array([Y]).T
#为自变量系数矩阵X赋值
X11=np.array([X1]).T
X22=np.array([X2]).T
A=np.array([[1],[1],[1],[1],[1],[1],[1],[1],[1],[1]])#创建系数矩阵
B=np.hstack((A,X11))#将矩阵a与矩阵X11合并为矩阵b
X=np.hstack((B,X22))#将矩阵b与矩阵X22合并为矩阵X
#求矩阵X的转置矩阵
X_=X.T
#求矩阵X与他的转置矩阵的X_的乘积
X_X=np.dot(X_,X)
#求矩阵X与他的转置矩阵的X_的乘积的逆矩阵
X_X_=np.linalg.inv(X_X)
#求解系数矩阵W，分别对应截距b、a1、和a2
W=np.dot(np.dot((X_X_),(X_)),Y1)
b=W[0][0]
a1=W[1][0]
a2=W[2][0]
print("系数a1={:.1f}".format(a1))
print("系数a2={:.1f}".format(a2))
print("截距为={:.1f}".format(b))
print("多元线性回归方程为:y={:.1f}".format(a1),"X1+ {:.1f}".format(a2),"X2+{:.1f}".format(b))
#画出线性回归分析图
data1=pd.read_excel('D:\面积距离车站数据.xlsx')
sns.pairplot(data1, x_vars=['area','distance'], y_vars='Y', height=3, aspect=0.8, kind='reg')
plt.show()
#求月销售量Y的和以及平均值y1
sumy=0#因变量的和
y1=0#因变量的平均值
for i in range(0,len(Y)):sumy=sumy+Y[i]
y1=sumy/len(Y)
#求月销售额y-他的平均值的和
y_y1=0#y-y1的值的和
for i in range(0,len(Y)):y_y1=y_y1+(Y[i]-y1)
print("营业额-营业额平均值的和为:{:.1f}".format(y_y1))
#求预测值sales1
sales1=[]
for i in range(0,len(Y)):sales1.append(a1*X1[i]+a2*X2[i]+b)
#求预测值的平均值y2
y2=0
sumy2=0
for i in range(len(sales1)):sumy2=sumy2+sales1[i]
y2=sumy2/len(sales1)
#求预测值-平均值的和y11_y2
y11_y2=0
for i in range(0,len(sales1)):y11_y2=y11_y2+(sales1[i]-y2)
print("预测营业额-预测营业额平均值的和为:{:.1f}".format(y11_y2))
#求月销售额y-他的平均值的平方和
Syy=0#y-y1的值的平方和
for i in range(0,len(Y)):Syy=Syy+((Y[i]-y1)*(Y[i]-y1))
print("Syy={:.1f}".format(Syy))
#求y1-y1平均的平方和
Sy1y1=0
for i in range(0,len(sales1)):Sy1y1=Sy1y1+((sales1[i]-y2)*(sales1[i]-y2))
print("Sy1y1={:.1f}".format(Sy1y1))
#（y1-y1平均）*（y-y平均）
Syy1=0
for i in range(0,len(sales1)):Syy1=Syy1+((Y[i]-y1)*(sales1[i]-y2))
print("Syy1={:.1f}".format(Syy1))
#求y-y1的平方Se
Se=0
for i in range(0,len(sales1)):Se=Se+((Y[i]-y2)*(Y[i]-y2))
print("Se={:.1f}".format(Se))
#求R
R=Syy1/((Syy*Sy1y1)**0.5)
R2=R*R
print("R2={:.4f}".format(R2))
运行结果
2、采用sklearn库
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from sklearn import linear_model
%matplotlib inlinedata=pd.read_excel('D:\面积距离车站数据.xlsx') #导入数据
X=data[['area','distance']] #x,y取值
y=data['Y']
y=data.Y
model = linear_model.LinearRegression()
model.fit(X,y)
print("系数a1={:.1f}".format(model.coef_[0]))
print("系数a2={:.1f}".format(model.coef_[1]))
print("截距为={:.1f}".format(model.intercept_))
print("多元线性回归方程为:y={:.1f}".format(model.coef_[0]),"X1+ {:.1f}".format(model.coef_[1]),"X2+{:.1f}".format(model.intercept_))
sns.pairplot(data, x_vars=['area','distance'], y_vars='Y', height=3, aspect=0.8, kind='reg')
plt.show()
运行结果
（四）将梯度下降法与最小二乘法进行对比
这里给出利用Excel表格的数据分析功能得出的数据：
通过对两种方法的对比分析，不难得出结论：最小二乘法计算量偏大，而且求逆矩阵相当耗费时间，求逆矩阵也会存在数值不稳定的情况，相比之下梯度下降法可以看做是一种更简单的最小二乘法最后一步解方程的方法，梯度下降法虽然有一定的缺点，但是计算量不大，在数据量比较大的时候选择梯度下降法比较好一点。
利用Python编程，分别使用梯度下降法和最小二乘法求解多元函数相关推荐
- 梯度下降法和Sklearn实现线性回归
1.线性回归模型 线性回归 (linear regression)是-种线性模型, 它假设输入变量 x和单个输出变量 y 之间存在线性关系.具体来说,利用线性回归模型,可以从-组输入变量 x 的线性组 ...
- python二元函数求导_用Excel和python实现二元函数梯度下降的人工智能,之用,excel,一元...
梯度下降法和牛顿法的总结与比较 机器学习的本质是建立优化模型,通过优化方法,不断迭代参数向量,找到使目标函数最优的参数向量.最终建立模型 通常用到的优化方法:梯度下降方法.牛顿法.拟牛顿法等.这些优化 ...
- 八、梯度下降法和拟牛顿法
1.梯度 2.梯度上升和梯度下降 3.梯度下降算法详解 3.1 直观解释 3.2 梯度下降相关概念 3.3 梯度下降的矩阵描述 3.4 梯度下降的算法调优 4.梯度下降法大家族 5.梯度下降法和其他无 ...
- 梯度下降算法_批梯度下降法，Minibatch梯度下降法和随机梯度下降法之间的区别...
什么是梯度下降法? 梯度下降法是一种机器学习中常用的优化算法,用来找到一个函数(f)的参数(系数)的值,使成本函数(cost)最小. 当参数不能解析计算时(如使用线性代数),并且必须通过优化算法搜索时 ...
- GBDT与xgb区别，以及梯度下降法和牛顿法的数学推导
为什么要介绍梯度下降法和牛顿法那? 这里提及两个算法模型GBDT和XGBoost,两个都是boosting模型. GBDT和xgb的目标函数是不同的,同时针对其目标函数中的误差函数 L(θ) 的拟合方 ...
- 【python】三种梯度下降学习率策略的比较(exact line search, backtracking, diminishing steps)
简要介绍: 1. exact line search 即在梯度下降的每次迭代中选择使梯度下降最大的学习率.我们可以使用黄金分割法来求解. 关于黄金分割法求根的实现golden_section()见专栏 ...
- Python：利用python编程实现三维图像绘制展示(六面体旋转、三维球柱状体、下雪场景等)
Python:利用python编程实现三维图像绘制展示(六面体旋转.三维球柱状体.下雪场景等) 目录 利用python编程实现三维图像绘制展示(六面体旋转.三维球柱状体.下雪场景等) 1.实现六面体旋 ...
- Python：利用python编程将上海十六区，2020年5月份房价实时地图(数据来源房天下)进行柱状图、热图可视化
Python:利用python编程将上海十六区,2020年5月份房价实时地图(数据来源房天下)进行柱状图.热图可视化 目录 上海十六区,2020年5月份房价实时地图(数据来源房天下)可视化 雷达图.柱 ...
- pyaudio：基于pyaudio利用Python编程从电脑端录制音频保存到指定文件夹+将录音上传服务器+录音进行识别并转为文本保存
pyaudio:基于pyaudio利用Python编程从电脑端录制音频保存到指定文件夹+将录音上传服务器+录音进行识别并转为文本保存 目录 输出结果 代码实现 输出结果 代码实现 # -*- codi ...
- 机器学习中梯度下降法和牛顿法的比较
在机器学习的优化问题中,梯度下降法和牛顿法是常用的两种凸函数求极值的方法,他们都是为了求得目标函数的近似解.在逻辑斯蒂回归模型的参数求解中,一般用改良的梯度下降法,也可以用牛顿法.由于两种方法有些相似 ...
最新文章
- Eclipse Color Theme
- oom 如何避免 高并发_糖尿病并发症真的会致死？又该如何避免它发生？
- 「Ubuntu」Ubuntu中的python终端配置（修改终端默认python配置，软连接，不同版本python环境配置）
- token 过期刷新令牌_OkHttp实现全局过期token自动刷新
- NSTimer不准确与GCDTimer详解
- android类似QQ空间，微信朋友圈，微博主页源码
- 阿里巴巴数据报告：消费向境内回流 低线城市消费蓬勃
- react学习系列1 修改create-react-app配置支持stylus 1
- 全栈溯源、mAPM、金融性能、Oracle VS. MySQL：看APM技术专场有哪些干货
- 初始化对于类与接口的异同点深入解析
- 关于某学习通网页鼠标不能移出视频窗口的问题
- 家用计算机 阵列,家用电脑如何建立RAID？
- 解决手机不能设置DeviceOwner权限提示already provisioned问题
- 分享一个盟重英雄脚本挂机工具（附随机数生成源码)
- 圣诞节用java画一棵圣诞树给你的女友
- 三菱FX3U/FX1N底层源码,PLSR源码， 总体功能和指令可能支持在RUN中下载程序，支持注释的写入和读取，支持脉冲输出与定位指令(包括PLSY /PWM/PLSR/PLSV/DRVI /DRV
- 架构师的 36 项修炼1 开篇词：7分钟Get技术人进阶技巧
- 4.1 用格雷戈里公式求π的近似值
- LoRa手持无线终端相比我们常用的工业PDA有哪些优势
- CSS 实现鼠标触碰完成渐变切换