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=ax^n+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

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