python最小二乘法拟合直线_python_numpy最小二乘法的曲线拟合

在了解了最小二乘法的基本原理之后python_numpy实用的最小二乘法理解，就可以用最小二乘法做曲线拟合了

1.直线拟合

直线拟合

已知图中拟合数据的坐标，对图中的拟合数据进行直线拟合。

依旧使用最小二乘法求解

Ax=b——————1

无解下的最优解。已知点的个数为n，所求直线的方程为y1=ax1+b，A由方程右边的a，b的系数构成构成(nx2)的矩阵,每行为(x1,1)，b由已知点的y1坐标构成矩阵(nx1)。方程1中的x为要求的列向量[a,b]。

A.TAx'=A.Tb

x'=(A.TA)^(-1)A.TC

求得x‘后，画出拟合曲线的yy=Ax'

import numpy as np

import matplotlib.pyplot as plt

#x的个数决定了样本量

x = np.arange(-1,1,0.02)

#y为理想函数

y = 2*np.sin(x*2.3)+0.5*x**3

#y1为离散的拟合数据

y1 = y+0.5*(np.random.rand(len(x))-0.5)

##################################

#主要程序

one=np.ones((len(x),1))#len(x)得到数据量

x=x.reshape(x.shape[0],1)

A=np.hstack((x,one))#两个100x1列向量合并成100x2,(100, 1) (100,1 ) (100, 2)

C=y1.reshape(y1.shape[0],1)

#等同于C=y1.reshape(100,1)

#虽然知道y1的个数为100但是程序中不应该出现人工读取的数据

def optimal(A,b):

B = A.T.dot(b)

AA = np.linalg.inv(A.T.dot(A))#求A.T.dot(A)的逆

P=AA.dot(B)

print P

return A.dot(P)

#求得的[a,b]=P=[[ 2.88778507e+00] [ -1.40062271e-04]]

yy = optimal(A,b)

#yy=P[0]*x+P[1]

##################################

plt.plot(x,y,color='g',linestyle='-',marker='',label=u'理想曲线')

plt.plot(x,y1,color='m',linestyle='',marker='o',label=u'拟合数据')

plt.plot(x,yy,color='b',linestyle='-',marker='.',label=u"拟合曲线")

# 把拟合的曲线在这里画出来

plt.legend(loc='upper left')

plt.show()

直线拟合结果

从结果中可以看出，直线拟合并不能对拟合数据达到很好的效果，下面我们介绍一下曲线拟合。

2.曲线拟合

曲线拟合

图中的拟合数据如果用直线进行拟合效果会更差，曲线能更好的表达数据的特征。这里我们使用多项式函数进行拟合。

拟合函数：

y=axn+bx(n-1)+cx^(n-2)+...+d

假设拟合数据共有100个

由 Ax=b

A=[x1^n x1^(n-1) x1^(n-2) ...... 1]

[x2^n x2^(n-1) x2^(n-2) ...... 1]

......

[x100^n x100^(n-1) x100^(n-2) . 1]

b=[y1]

[y2]

......

[y100]

解得拟合函数的系数[a,b,c.....d]

CODE:

import numpy as np

import matplotlib.pyplot as plt

x = np.arange(-1,1,0.02)

y = ((x*x-1)**3+1)*(np.cos(x*2)+0.6*np.sin(x*1.3))

y1 = y+(np.random.rand(len(x))-0.5)

##################################

### 核心程序

#使用函数y=ax^3+bx^2+cx+d对离散点进行拟合，最高次方需要便于修改，所以不能全部列举，需要使用循环

#A矩阵

m=[]

for i in xrange(7):#这里选的最高次为x^7的多项式

a=x**(i)

m.append(a)

A=np.array(m).T

b=y1.reshape(y1.shape[0],1)

##################################

def projection(A,b):

AA = A.T.dot(A)#A乘以A转置

w=np.linalg.inv(AA).dot(A.T).dot(b)

print w#w=[[-0.03027851][ 0.1995869 ] [ 2.43887827] [ 1.28426472][-5.60888682] [-0.98754851][ 2.78427031]]

return A.dot(w)

yw = projection(A,b)

yw.shape = (yw.shape[0],)

plt.plot(x,y,color='g',linestyle='-',marker='',label=u"理想曲线")

plt.plot(x,y1,color='m',linestyle='',marker='o',label=u"已知数据点")

plt.plot(x,yw,color='r',linestyle='',marker='.',label=u"拟合曲线")

plt.legend(loc='upper left')

plt.show()

结果

根据结果可以看到拟合的效果不错。

我们可以通过改变

拟合函数类型

样本数(此处为x的个数)

来调整拟合效果。

如果此处我们把拟合函数改为最高次为x^20的多项式

m=[]

for i in xrange(20):

a=x**(i)

m.append(a)

所得结果如下：

x^20 样本数100

这种现象称为过拟合现象

可以通过增加样本数数，

降低拟合函数的次数

矫正过拟合现象

在保持拟合函数改为最高次为x^20的多项式的条件下，增大样本数：

x = np.arange(-1,1,0.005) #原来是x = np.arange(-1,1,0.02)

x^20 样本数400

通过结果可以看出，过拟合现象得到了改善。