用Excel和Python编程，实现梯度法求解二次函数极小值和漫画书中“店铺多元回归”问题（还有牛顿法）

用Excel和Python编程，实现梯度法求解二次函数极小值和漫画书中“店铺多元回归”问题（还有牛顿法）

梯度下降法定义：

梯度下降算法最速下降法又称为梯度法，它是解析法中最古老的一种，其他解析方法或是它的变形，或是受它的启发而得到的，因此它是最优化方法的基础。作为一种基本的算法，他在最优化方法中占有重要地位

牛顿法：

在最优化的问题中，线性最优化至少可以使用单纯行法求解，但对于非线性优化问题，牛顿法提供了一种求解的办法。假设任务是优化一个目标函数f，求函数f的极大极小问题，可以转化为求解函数f的导数f’=0的问题，这样求可以把优化问题看成方程求解问题（f’=0）。剩下的问题就和第一部分提到的牛顿法求解很相似了。

这次为了求解f’=0的根，把f（x）的泰勒展开，展开到2阶形式：

这个式子是成立的，当且仅当 Δx 无线趋近于0。此时上式等价与：

求解：

得出迭代公式：

一般认为牛顿法可以利用到曲线本身的信息，比梯度下降法更容易收敛（迭代更少次数）

二元函数的求解 :

二次函数$$f(x_1,x_2)=x_1^2+2x_2^2-4x_1-2x_1{x}_2$$
初始点$$[1,1{]}^T$$
迭代精度 $$\mathfrak{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+2x_2^2-4x_1-2x_1{x}_2$$
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.youkuaiyun.com/qq_35946969/article/details/84446000
[4]: http://www.wangyang0821.club/?p=90