python二维数据点的曲线拟合

这里采用最小二乘法对数据进行拟合，x的最高幂次可调

即：

原理见隔壁大佬的“最小二乘法 来龙去脉” 

代码：

import matplotlib.pylab as pl
from numpy import array, linalg
import math as m
import numpy as np

 

读取文件

#---------------------------------------------读取文件
xlist = []
ylist = []

istep=0
with open ('cal4.txt', 'r') as infile:
     for lines in infile:
         words=lines.split()
         xlist.append(float(words[0]))
         ylist.append(float(words[1]))
         istep = istep + 1
nstep = istep

xy坐标储存在xlist,ylist中，nstep是坐标数量

输入最高次数n

#------------------------------------------------输入最高次数
line = input("input the maximal x^n:")
n = int(line)
print("the maximal number of n is:", n)
print("the maximal number of walk steps is:", nstep)

 

求矩阵A

见上述链接，为求系数矩阵β，需求数据x的次幂矩阵A，

#------------------------------------------------求系数矩阵
#初始化矩阵A
x0 = xlist[0]
xarray0 = []
for i in range(0,n+1):
    xarray0.append(m.pow(x0,i))
A = array([xarray0])

#求矩阵A
for xn in range(1,nstep):              #因为在初始化中已经把X1记入A，因此这里从xlist[1]即X2开始
    x = xlist[xn]
    arx = []
    for i in range(0,n+1):             #对于每个Xm，得到其次幂列表[1，X，X^2，...，X^n]
        arx.append(m.pow(x,i))         
 
    arrx = array([arx])                    #把列表arx记为一维矩阵arrx，否则会报错
    A = np.concatenate((A,arrx),axis=0)    #通过concatenate函数将arrx合入A的行

如果直接令A = array([])，concatenate会报错：维数不符合。所以找到的一个繁琐的解决办法是先求出X1的次幂列表，用它初始化A，后续的维数就能够相符合了

用A和Y求系数矩阵β

AT = A.T
B = np.dot(AT,A)
C = linalg.inv(B)
D = np.dot(C,AT)

Y = []
for y in ylist:
    Y.append(y)

beita = np.dot(D,Y)

print("其系数由低到高为：", beita)

画图

记得按照取点范围调整x的画图范围      xest = np.arange(min, max, 0.1)

#-----------------------------------------------------------画图
#列出函数式 y = y(x)
xest = np.arange(0, 15, 0.1)       #这里的x范围需要自己调整，也可以再添加程序令其覆盖xmin到xmax
yest = 0
xplus = 1

for i in range(0,n+1):
    yest = yest + beita[i] * xplus
    xplus *= xest


#画曲线
pl.xlabel('x')
pl.ylabel('y')
pl.title("OK")
pl.plot(xest, yest)

#画分布的点
pl.plot(xlist, ylist, 'ro', lw=2, markersize=6)
pl.show()

效果：

初学python，有很多地方只会用笨办法避免报错，有更好的地方欢迎指正