數值分析實驗之最小二乘擬合 含有噪聲擾動（python實現）

一、實驗目的

  掌握最小二乘法擬合離散數據，多項式函數形式擬合曲線以及可以其他可以通過變量變換轉化為多項式的擬合曲線目前待實現功能：

   1. 最小二乘法的基本實現。

   2. 用不同數據量，不同參數，不同的多項式階數，比較實驗效果。

   3. 語言python。

二、實驗原理

  最小二乘法（又稱最小平方法）是一種數學優化技術。它通過最小化誤差的平方和尋找數據的最佳函數匹配。利用最小二乘法可以簡便地求得未知的數據，並使得這些求得的數據與實際數據之間誤差的平方和為最小。最小二乘法還可用於曲線擬合。其他一些優化問題也可通過最小化能量或最大化熵用最小二乘法來表達。

三、實驗內容

  求y=f(x)=sin(x)+h(x)在區間[0,10]上按101等距節點確定的離散數據點組(x_i,y_i)的直線擬合以及曲線擬合，其中是服從h(x)標准正態分布的噪聲擾動

四、程序實現

 • 一次擬合：

import numpy as np
import matplotlib.pyplot as plt
import math
#定義x、y散點坐標
x = np.arange(0.0, 10.0,0.1)
x = np.array(x)
print('x is :\n',x)
num = np.sin(x)+np.random.randn(100)
y = np.array(num)
print('y is :\n',y)
f1 = np.polyfit(x, y, 1)#用1次多項式擬合，若要多次擬合，相應的改變這個常數即可
print('f1 is :\n',f1)

p1 = np.poly1d(f1)
print('p1 is :\n',p1)

#也可使用yvals=np.polyval(f1, x)
yvals = p1(x)  #擬合y值
print('yvals is :\n',yvals)
#繪圖
plot1 = plt.plot(x, y, 's',label='original values',color="blue")
plot2 = plt.plot(x, yvals, 'r',label='polyfit values',color="red")
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('polyfitting')
plt.show()

   運行結果：

              

   所得圖形：

   

   • 曲線擬合（用a*sin(x)+b擬合）：

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy.optimize import curve_fit

#自定義函數 
def func(x, a, b):
    return a*np.sin(x)+b

#定義x、y散點坐標
x = np.arange(0.0, 10.0,0.1)
x = np.array(x)
num = np.sin(x)+np.random.randn(100)
y = np.array(num)

#非線性最小二乘法擬合
popt, pcov = curve_fit(func, x, y)
#獲取popt里面是擬合系數
print(popt)
a = popt[0]
b = popt[1]
#c = popt[2]
#d = popt[3]
#e = popt[4]
yvals = func(x,a,b) #擬合y值
print('popt:', popt)
print('系數a:', a)
print('系數b:', b)
#print('系數c:', c)
#print('系數d:', d)
#print('系數e:', e)
print('系數pcov:', pcov)#方差
print('系數yvals:', yvals)#x代入擬合出的函數得到的函數值
#繪圖
plot1 = plt.plot(x, y, 's',label='original values',color="purple")
x_test = np.arange(0.0, 10.0, 0.01)
y_test = func(x_test,a,b)
plot2 = plt.plot(x_test, y_test, 'r',label='polyfit values',color="red")
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

運行結果

                

   所得圖形：

     

  •曲線擬合（用a*np.sin(b*x+c)+d擬合）：

 

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy.optimize import curve_fit

#自定義函數 
def func(x, a, b, c, d):
    return a*np.sin(b*x+c)+d

#定義x、y散點坐標
x = np.arange(0.0, 10.0,0.1)
x = np.array(x)
num = np.sin(x)+np.random.randn(100)
y = np.array(num)

#非線性最小二乘法擬合
popt, pcov = curve_fit(func, x, y)
#獲取popt里面是擬合系數
print(popt)
a = popt[0]
b = popt[1]
c = popt[2]
d = popt[3]
yvals = func(x,a,b,c,d) #擬合y值
print('popt:', popt)
print('系數a:', a)
print('系數b:', b)
print('系數c:', c)
print('系數d:', d)
print('系數pcov:', pcov)#方差
print('系數yvals:', yvals)#x代入擬合出的函數得到的函數值
#繪圖
plot1 = plt.plot(x, y, 's',label='original values',color='orange')
x_test = np.arange(0.0, 10.0, 0.01)
y_test = func(x_test,a,b,c,d)
plot2 = plt.plot(x_test, y_test, 'r',label='polyfit values',color='brown')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

   運行結果：

              

   所得圖形：

   

  •自定義函數實現：

 

import numpy as np
import matplotlib.pyplot as plt
import math
from scipy.optimize import curve_fit

#自定義函數 
def func(x, a0,a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8,w):
    return  a0 + a1*np.cos(x*w) + b1*np.sin(x*w) + \
               a2*np.cos(2*x*w) + b2*np.sin(2*x*w) + a3*np.cos(3*x*w) + b3*np.sin(3*x*w) + \
               a4*np.cos(4*x*w) + b4*np.sin(4*x*w) + a5*np.cos(5*x*w) + b5*np.sin(5*x*w) + \
               a6*np.cos(6*x*w) + b6*np.sin(6*x*w) + a7*np.cos(7*x*w) + b7*np.sin(7*x*w) + \
               a8*np.cos(8*x*w) + b8*np.sin(8*x*w)

#定義x、y散點坐標
x = np.arange(0.0, 10.0,0.1)
x = np.array(x)
num = np.sin(x)+np.random.randn(100)
y = np.array(num)

#非線性最小二乘法擬合
popt, pcov = curve_fit(func, x, y)
#獲取popt里面是擬合系數
print(popt)
a0 = popt[0]
a1 = popt[1]
a2 = popt[2]
a3 = popt[3]
a4 = popt[4]
a5 = popt[5]
a6 = popt[6]
a7 = popt[7]
a8 = popt[8]
b1 = popt[9]
b2 = popt[10]
b3 = popt[11]
b4 = popt[12]
b5 = popt[13]
b6 = popt[14]
b7 = popt[15]
b8 = popt[16]
w = popt[17]
yvals = func(x,a0,a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8,w) #擬合y值
print('popt:', popt)
print('系數a0:', a0)
print('系數a1:', a1)
print('系數a2:', a2)
print('系數a3:', a3)
print('系數a4:', a4)
print('系數a5:', a5)
print('系數a6:', a6)
print('系數a7:', a7)
print('系數a8:', a8)
print('系數b1:', b1)
print('系數b2:', b2)
print('系數b3:', b3)
print('系數b4:', b4)
print('系數b5:', b5)
print('系數b6:', b6)
print('系數b7:', b7)
print('系數b8:', b8)
print('系數w:', w)
print('系數pcov:', pcov)#方差
print('系數yvals:', yvals)#x代入擬合出的函數得到的函數值
#繪圖
plot1 = plt.plot(x, y, 's',label='original values',color='yellow')
x_test = np.arange(0.0, 10.0, 0.01)
y_test = func(x_test,a0,a1,a2,a3,a4,a5,a6,a7,a8,b1,b2,b3,b4,b5,b6,b7,b8,w)
plot2 = plt.plot(x_test, y_test, 'r',label='polyfit values',color='blue')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()

   實現結果：

           

   所得圖形：

    

心得體會

  通過本次實驗，我對MATLAB的操作更加熟悉，也對本學期正在學習的Python有了更深層次的認識，對着兩種編程軟件更加熟悉。最重要的是：對最小二乘法理解深入，以及可以應用推廣，對數值分析理論的學習有很重要的作用，總而言之收獲頗多。