在了解了最小二乘法的基本原理之后python_numpy實用的最小二乘法理解饶深,就可以用最小二乘法做曲線擬合了
1.直線擬合
直線擬合
已知圖中擬合數據的坐標数尿,對圖中的擬合數據進行直線擬合蜜另。
依舊使用最小二乘法求解
Ax=b——————1
無解下的最優(yōu)解。已知點的個數為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對離散點進行擬合,最高次方需要便于修改节值,所以不能全部列舉徙硅,需要使用循環(huán)
#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
通過結果可以看出匿乃,過擬合現象得到了改善桩皿。