用Excel和Python编程,实现梯度法求解二次函数极小值和漫画书中“店铺多元回归”问题(还有牛顿法)
梯度下降法定义:
梯度下降算法最速下降法又称为梯度法,它是解析法中最古老的一种,其他解析方法或是它的变形,或是受它的启发而得到的,因此它是最优化方法的基础。作为一种基本的算法,他在最优化方法中占有重要地位
牛顿法:
在最优化的问题中,线性最优化至少可以使用单纯行法求解,但对于非线性优化问题,牛顿法提供了一种求解的办法。假设任务是优化一个目标函数f,求函数f的极大极小问题,可以转化为求解函数f的导数f’=0的问题,这样求可以把优化问题看成方程求解问题(f’=0)。剩下的问题就和第一部分提到的牛顿法求解很相似了。
这次为了求解f’=0的根,把f(x)的泰勒展开,展开到2阶形式:
这个式子是成立的,当且仅当 Δx 无线趋近于0。此时上式等价与:
求解:
得出迭代公式:
一般认为牛顿法可以利用到曲线本身的信息,比梯度下降法更容易收敛(迭代更少次数)
二元函数的求解 :
二次函数
f
(
x
1
,
x
2
)
=
x
1
2
+
2
x
2
2
−
4
x
1
−
2
x
1
x
2
f(x_1,x_2)=x_1^2+2x_2^2-4x_1-2x_1x_2
f(x1,x2)=x12+2x22−4x1−2x1x2
初始点
[
1
,
1
]
T
[1,1]^T
[1,1]T
迭代精度
y
=
0.001
\mathfrak y=0.001
y=0.001
1.Excel:
分别求出x1,x2的偏导是关键
根据算法求出极值点为-7.99999
2.Python编程:
要先导入包
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
from mpl_toolkits.mplot3d import Axes3D
import warnings
防止中文乱码
import matplotlib as mpl
# 设置字符集,防止中文乱码
mpl.rcParams['font.sans-serif'] = [u'simHei']
mpl.rcParams['axes.unicode_minus'] = False
1、构建一个二次函数
f
(
x
1
,
x
2
)
=
x
1
2
+
2
x
2
2
−
4
x
1
−
2
x
1
x
2
f(x_1,x_2)=x_1^2+2x_2^2-4x_1-2x_1x_2
f(x1,x2)=x12+2x22−4x1−2x1x2
2、随机生成X1,X2点,根据X1,X2点生成Y点。
3、画出图像。
def f2(x1,x2):
return x1**2+2*x2**2-4*x1-2*x1*x2
X1 = np.arange(-4,4,0.2)
X2 = np.arange(-4,4,0.2)
X1, X2 = np.meshgrid(X1, X2) # 生成xv、yv,将X1、X2变成n*m的矩阵,方便后面绘图
Y = np.array(list(map(lambda t : f2(t[0],t[1]),zip(X1.flatten(),X2.flatten()))))
Y.shape = X1.shape # 1600的Y图还原成原来的(40,40)
%matplotlib inline
#作图
fig = plt.figure(facecolor='w')
ax = Axes3D(fig)
ax.plot_surface(X1,X2,Y,rstride=1,cstride=1,cmap=plt.cm.jet)
ax.set_title(u'$ y = x1^2+2*x2^2-4*x1-2*x1*x2 $')
plt.show()
图如下:
梯度法:
对当前二维图像求最小点
1、随机取一个点(x1,x2),设定α参数值
2、对这个点分别求关于x1、x2的偏导数
# 二维原始图像
def f2(x, y):
return x**2+2*(y)**2 - 4*(x)-2*(x)*(y)
## 偏函数
def hx1(x, y):
return 2*x-4-2*y
def hx2(x, y):
return 4*y-2*x
X1 = np.arange(-4,4,0.2)
X2 = np.arange(-4,4,0.2)
X1, X2 = np.meshgrid(X1, X2) # 生成xv、yv,将X1、X2变成n*m的矩阵,方便后面绘图
Y = np.array(list(map(lambda t : f2(t[0],t[1]),zip(X1.flatten(),X2.flatten()))))
Y.shape = X1.shape # 1600的Y图还原成原来的(40,40)
x1 = 1
x2 = 1
alpha = 0.1
#保存梯度下降经过的点
GD_X1 = [x1]
GD_X2 = [x2]
GD_Y = [f2(x1,x2)]
# 定义y的变化量和迭代次数
y_change = f2(x1,x2)
iter_num = 0
while(y_change < 1e-10 and iter_num < 100) :
tmp_x1 = x1 - alpha * hx1(x1,x2)
tmp_x2 = x2 - alpha * hx2(x1,x2)
tmp_y = f2(tmp_x1,tmp_x2)
f_change = np.absolute(tmp_y - f2(x1,x2))
x1 = tmp_x1
x2 = tmp_x2
GD_X1.append(x1)
GD_X2.append(x2)
GD_Y.append(tmp_y)
iter_num += 1
print(u"最终结果为:(%.5f, %.5f, %.5f)" % (x1, x2, f2(x1,x2)))
print(u"迭代过程中X的取值,迭代次数:%d" % iter_num)
print(GD_X1)
再作图:
# 作图
fig = plt.figure(facecolor='w',figsize=(20,18))
ax = Axes3D(fig)
ax.plot_surface(X1,X2,Y,rstride=1,cstride=1,cmap=plt.cm.jet)
ax.plot(GD_X1,GD_X2,GD_Y,'ko-')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.set_title(u'函数;\n学习率:%.3f; 最终解:(%.3f, %.3f, %.3f);迭代次数:%d' % (alpha, x1, x2, f2(x1,x2), iter_num))
plt.show()
漫画书中“店铺多元回归”问题
1.梯度法:
导入:
import numpy as np
from matplotlib import pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
data=np.genfromtxt('D:\\重回归.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=0
for i in range(0,len(x_data)):#定义一共有多少样本点
totalerror=totalerror+(y_data[i]-(theta1*x_data[i,0]+theta2*x_data[i,1]+theta0))**2
return 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=0
theta1_grad=0
theta2_grad=0
for 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_grad
theta1=theta1-lr*theta1_grad
theta2=theta2-lr*theta2_grad
return 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()
结果如图:
2.牛顿法:
1.导入库
from sklearn import linear_model
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
导入数据,并求出x1,x2,d
data=pd.read_excel('D:\\重回归.xlsx')
X=data[['X1','X2']]
Y=data['Y']
Y=data.Y
model = linear_model.LinearRegression()
model.fit(X,Y)
a1=model.coef_[0]
a2=model.coef_[1]
b=model.intercept_
print("X1=",a1)
print("X2=",a2)
print("d=",b)
print("多元线性回归方程为:y=",a1,"X1",a2,"X2+",b)
sns.pairplot(data, x_vars=['X1','X2'], 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("y-y的平均数的和为:",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("(^y-^y的平均数)的和为:",y11_y2)
#求月销售额y-y的平均数的平方和
Syy=0#y-y1的值的平方和
for i in range(0,len(Y)):
Syy=Syy+((Y[i]-y1)*(Y[i]-y1))
print("Syy=",Syy)
#求y1-y1平均的平方和即(^y-^y的平均数)的平方
Sy1y1=0
for i in range(0,len(sales1)):
Sy1y1=Sy1y1+((sales1[i]-y2)*(sales1[i]-y2))
print("Sy1y1=",Sy1y1)
#求(y1-y1平均)*(y-y的平均数)
Syy1=0
for i in range(0,len(sales1)):
Syy1=Syy1+((Y[i]-y1)*(sales1[i]-y2))
print("Syy1=",Syy1)
#求R
R=Syy1/((Syy*Sy1y1)**0.5)
R2=R*R
print("R=",R)
print("判定系数R^2=",R2)
参考文献
[1].基于jupyter notebook的python编程-----通过矩阵实现最小二乘法对求多元线性回归方程的的待定系数a和判定系数R2的求解
[2]线性回归 - 多元线性回归案例 - 分析步骤、输出结果详解、与Python的结果对比 -(SPSS建模)
[3]: https://blog.csdn.net/qq_35946969/article/details/84446000
[4]: http://www.wangyang0821.club/?p=90