對於樣本數據的散點圖形如函數y=ax2+bx+c的圖像的數據, 在python中的擬合過程為:
##最小二乘法 import numpy as np import scipy as sp import matplotlib.pyplot as plt from scipy.optimize import leastsq ''' 設置樣本數據,真實數據需要在這里處理 ''' ##樣本數據(Xi,Yi),需要轉換成數組(列表)形式 Xi=np.array([1,2,3,4,5,6]) #Yi=np.array([9,18,31,48,69,94]) Yi=np.array([9.1,18.3,32,47,69.5,94.8]) ''' 設定擬合函數和偏差函數 函數的形狀確定過程: 1.先畫樣本圖像 2.根據樣本圖像大致形狀確定函數形式(直線、拋物線、正弦余弦等) ''' ##需要擬合的函數func :指定函數的形狀 def func(p,x): a,b,c=p return a*x*x+b*x+c ##偏差函數:x,y都是列表:這里的x,y更上面的Xi,Yi中是一一對應的 def error(p,x,y): return func(p,x)-y ''' 主要部分:附帶部分說明 1.leastsq函數的返回值tuple,第一個元素是求解結果,第二個是求解的代價值(個人理解) 2.官網的原話(第二個值):Value of the cost function at the solution 3.實例:Para=>(array([ 0.61349535, 1.79409255]), 3) 4.返回值元組中第一個值的數量跟需要求解的參數的數量一致 ''' #k,b的初始值,可以任意設定,經過幾次試驗,發現p0的值會影響cost的值:Para[1] p0=[10,10,10] #把error函數中除了p0以外的參數打包到args中(使用要求) Para=leastsq(error,p0,args=(Xi,Yi)) #讀取結果 a,b,c=Para[0] print("a=",a,"b=",b,"c=",c) print("cost:"+str(Para[1])) print("求解的擬合直線為:") print("y="+str(round(a,2))+"x*x+"+str(round(b,2))+"x+"+str(c)) ''' 繪圖,看擬合效果. matplotlib默認不支持中文,label設置中文的話需要另行設置 如果報錯,改成英文就可以 ''' #畫樣本點 plt.figure(figsize=(8,6)) ##指定圖像比例: 8:6 plt.scatter(Xi,Yi,color="green",label="樣本數據",linewidth=2) #畫擬合直線 x=np.linspace(0,12,100) ##在0-15直接畫100個連續點 y=a*x*x+b*x+c ##函數式 plt.plot(x,y,color="red",label="擬合直線",linewidth=2) plt.legend() #繪制圖例 plt.show()
運行結果:
a= 2.06607141425 b= 2.5975001036 c= 4.68999985496 cost:1 求解的擬合直線為: y=2.07x*x+2.6x+4.68999985496
在R中的擬合過程:(在控制台直接敲入或者放入腳本都可以)
###設置函數形式
func<-function(a,b,c){
a*x*x+b*x+c
}
###設置樣本數據
x<-c(1,2,3,4,5,6)
y<-c(9.1,18.3,32,47,69.5,94.8)
###把樣本數據轉換為符合nls函數需要的格式
d<-data.frame(y,x)
###執行求解過程:如果x,y值完全一一對應,匯報錯誤(循環次數超過了50這個最大值)
nlmod<-nls(y ~ func(a1,b1,c1),data=d,start=list(a1=1,b1=1,c1=1),trace=F)
###分析結果
summary(nlmod)
運行結果: