python最小二乘法 实现 曲面拟合

问题：

给定若干个三维离散点，任意求出一个方程使得尽可能的将这些点包含在其中，并用python画出函数图像

思路：

最小二乘法，根据最优性表示出一个表达式E，为了让E最小，那么对任意一个系数求偏导后的值应该为0，列出若干条式子后，再用内置的高斯消元，求解系数。百度文库中有若干篇论文描述最小二乘法的实现方法。

代码：

%matplotlib inline
from matplotlib import pyplot as plt
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

def fun(x): #处理符号问题
    round(x,2)
    if x>=0: return '+'+str(x)
    else: return str(x)

#主函数
if __name__ == '__main__':
    n=int(input("请输入一个n，代表点的个数"))
    #print(n)
    X=[]
    Y=[]
    Z=[]
    #读入n个点的坐标
    for i in range(n):
        a=[int(i) for i in input('please input: ').split()]
        X.append(a[0])
        Y.append(a[1])
        Z.append(a[2])
    
    #print(X)    
    #print(Y)
    #print(Z)
    
   #fig = plt.figure() #建立一个新的图像
   #ax = Axes3D(fig) #建立一个三维坐标系
   #ax.scatter(X,Y,Z,c='b')
   #ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap='rainbow')
   # plt.show()


    #求方程系数
    sigma_x = 0
    for i in X : sigma_x+=i
    sigma_y = 0
    for i in Y : sigma_y+=i
    sigma_z = 0
    for i in Z : sigma_z+=i
    sigma_x2 = 0
    for i in X : sigma_x2+=i*i
    sigma_y2 = 0
    for i in Y : sigma_y2+=i*i
    sigma_x3 = 0
    for i in X : sigma_x3+=i*i*i
    sigma_y3 = 0
    for i in Y : sigma_y3+=i*i*i
    sigma_x4 = 0
    for i in X : sigma_x4+=i*i*i*i
    sigma_y4 = 0
    for i in Y : sigma_y4+=i*i*i*i
    sigma_x_y = 0
    for i in range(n) : 
        sigma_x_y+=X[i]*Y[i]
    #print(sigma_xy)
    sigma_x_y2 = 0
    for i in range(n) : sigma_x_y2+=X[i]*Y[i]*Y[i]
    sigma_x_y3 = 0
    for i in range(n) : sigma_x_y3+=X[i]*Y[i]*Y[i]*Y[i]
    sigma_x2_y = 0
    for i in range(n) : sigma_x2_y+=X[i]*X[i]*Y[i]
    sigma_x2_y2 = 0
    for i in range(n) : sigma_x2_y2+=X[i]*X[i]*Y[i]*Y[i]
    sigma_x3_y = 0
    for i in range(n) : sigma_x3_y+=X[i]*X[i]*X[i]*Y[i]
    sigma_z_x2 = 0
    for i in range(n) : sigma_z_x2+=Z[i]*X[i]*X[i]
    sigma_z_y2 = 0
    for i in range(n) : sigma_z_y2+=Z[i]*Y[i]*Y[i]
    sigma_z_x_y = 0
    for i in range(n) : sigma_z_x_y+=Z[i]*X[i]*Y[i]
    sigma_z_x = 0
    for i in range(n) : sigma_z_x+=Z[i]*X[i]
    sigma_z_y = 0
    for i in range(n) : sigma_z_y+=Z[i]*Y[i]
    #print("-----------------------")
    #给出对应方程的矩阵形式
    a=np.array([[sigma_x4,sigma_x3_y,sigma_x2_y2,sigma_x3,sigma_x2_y,sigma_x2],
               [sigma_x3_y,sigma_x2_y2,sigma_x_y3,sigma_x2_y,sigma_x_y2,sigma_x_y],
               [sigma_x2_y2,sigma_x_y3,sigma_y4,sigma_x_y2,sigma_y3,sigma_y2],
               [sigma_x3,sigma_x2_y,sigma_x_y2,sigma_x2,sigma_x_y,sigma_x],
               [sigma_x2_y,sigma_x_y2,sigma_y3,sigma_x_y,sigma_y2,sigma_y],
               [sigma_x2,sigma_x_y,sigma_y2,sigma_x,sigma_y,n]])
    b=np.array([sigma_z_x2,sigma_z_x_y,sigma_z_y2,sigma_z_x,sigma_z_y,sigma_z])
    #高斯消元解线性方程
    res= np.linalg.solve(a,b)
    #print(a)
    #print(b)
    #print(x)
   # print("-----------------------")

    #输出方程形式
    print("z=%.6s*x^2%.6s*xy%.6s*y^2%.6s*x%.6s*y%.6s"%(fun(res[0]),fun(res[1]),fun(res[2]),fun(res[3]),fun(res[4]),fun(res[5])))
    #画曲面图和离散点
    fig = plt.figure()#建立一个空间
    ax =fig.add_subplot(111,projection='3d')# 3D坐标

    n = 256
    u = np.linspace(-20,20,n) # 创建一个等差数列
    x,y = np.meshgrid(u,u) #转化成矩阵

    #给出方程
    z=res[0]*x*x+res[1]*x*y+res[2]*y*y+res[3]*x+res[4]*y+res[5]
    #画出曲面
    ax.plot_surface(x,y,z,rstride=3,cstride=3,cmap=cm.jet)
    #画出点
    ax.scatter(X, Y, Z,c='r')
    

参考网站：

https://blog.csdn.net/qq_39735236/article/details/79066032

https://www.cnblogs.com/NanShan2016/p/5493429.html

http://www.cnblogs.com/jukan/p/7826851.html

https://blog.csdn.net/qq5q13638/article/details/78398997

https://blog.csdn.net/eddy_zheng/article/details/48713449