【python資料探勘課程】十四.Scipy呼叫curve_fit實現曲線擬合

        前面系列文章講過各種知識，包括繪製曲線、散點圖、冪分佈等，而如何在在散點圖一堆點中擬合一條直線，也變得非常重要。這篇文章主要講述呼叫Scipy擴充套件包的curve_fit函式實現曲線擬合，同時計算出擬合的函式、引數等。希望文章對你有所幫助，如果文章中存在錯誤或不足之處，還請海涵~

        前文推薦：
       【Python資料探勘課程】一.安裝Python及爬蟲入門介紹
       【Python資料探勘課程】二.Kmeans聚類資料分析及Anaconda介紹
       【Python資料探勘課程】三.Kmeans聚類程式碼實現、作業及優化
       【Python資料探勘課程】四.決策樹DTC資料分析及鳶尾資料集分析
       【Python資料探勘課程】五.線性迴歸知識及預測糖尿病例項
       【Python資料探勘課程】六.Numpy、Pandas和Matplotlib包基礎知識
       【Python資料探勘課程】七.PCA降維操作及subplot子圖繪製
       【Python資料探勘課程】八.關聯規則挖掘及Apriori實現購物推薦
       【Python資料探勘課程】九.迴歸模型LinearRegression簡單分析氧化物資料
       【python資料探勘課程】十.Pandas、Matplotlib、PCA繪圖實用程式碼補充
       【python資料探勘課程】十一.Pandas、Matplotlib結合SQL語句視覺化分析
       【python資料探勘課程】十二.Pandas、Matplotlib結合SQL語句對比圖分析
       【python資料探勘課程】十三.WordCloud詞雲配置過程及詞頻分析

一. Scipy介紹

        SciPy (pronounced "Sigh Pie") 是一個開源的數學、科學和工程計算包。它是一款方便、易於使用、專為科學和工程設計的Python工具包，包括統計、優化、整合、線性代數模組、傅立葉變換、訊號和影象處理、常微分方程求解器等等。
        官方地址：https://www.scipy.org/

        Scipy常用的模組及功能如下圖所示：
        強烈推薦劉神的文章：Scipy高階科學計算 - 劉一痕

二. 曲線擬合

1.多項式擬合

        首先通過numpy.arange定義x、y座標，然後呼叫polyfit()函式進行3次多項式擬合，最後呼叫Matplotlib函式進行散點圖繪製(x,y)座標，並繪製預測的曲線。
        完整程式碼：

#encoding=utf-8  
import numpy as np
import matplotlib.pyplot as plt

#定義x、y散點座標
x = np.arange(1, 16, 1)
num = [4.00, 5.20, 5.900, 6.80, 7.34,
       8.57, 9.86, 10.12, 12.56, 14.32,
       15.42, 16.50, 18.92, 19.58, 20.00]
y = np.array(num)

#用3次多項式擬合
f1 = np.polyfit(x, y, 3)
p1 = np.poly1d(f1)
print(p1)

#也可使用yvals=np.polyval(f1, x)
yvals = p1(x)  #擬合y值

#繪圖
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('polyfitting')
plt.show()
plt.savefig('test.png')

2.e的b/x次方擬合

        下面採用Scipy的curve_fit()對上面的資料進行e的b/x次方擬合。資料集如下：

x = np.arange(1, 16, 1)
num = [4.00, 5.20, 5.900, 6.80, 7.34,
       8.57, 9.86, 10.12, 12.56, 14.32,
       15.42, 16.50, 18.92, 19.58, 20.00]
y = np.array(num)

#encoding=utf-8  
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

#自定義函式 e指數形式
def func(x, a, b):
    return a*np.exp(b/x)

#定義x、y散點座標
x = np.arange(1, 16, 1)
num = [4.00, 5.20, 5.900, 6.80, 7.34,
       8.57, 9.86, 10.12, 12.56, 14.32,
       15.42, 16.50, 18.92, 19.58, 20.00]
y = np.array(num)

#非線性最小二乘法擬合
popt, pcov = curve_fit(func, x, y)
#獲取popt裡面是擬合係數
a = popt[0] 
b = popt[1]
yvals = func(x,a,b) #擬合y值
print u'係數a:', a
print u'係數b:', b

#繪圖
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.show()
plt.savefig('test2.png')

3.aX的b次方擬合

        第三種方法是通過Pandas匯入資料，因為通常資料都會儲存在csv、excel或資料庫中，所以這裡結合讀寫資料繪製a*x的b次方形式。
        假設本地存在一個data.csv檔案，資料集如下圖所示：

#匯入資料及x、y散點座標
data = pd.read_csv("data.csv")
print data
print(data.shape)    
print(data.head(5)) #顯示前5行資料
x = data['x'] #獲取x列
y = data['y'] #獲取y列
print x
print y

0      4.00
1      5.20
2      5.90
3      6.80
4      7.34
5      8.57
6      9.86
7     10.12
8     12.56
9     14.32
10    15.42
11    16.50
12    18.92
13    19.58
14    20.00
Name: y, dtype: float64

#encoding=utf-8  
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
import pandas as pd  

#自定義函式 e指數形式
def func(x, a, b):
    return a*pow(x,b)

#匯入資料及x、y散點座標
data = pd.read_csv("data.csv")
print data
print(data.shape)    
print(data.head(5)) #顯示前5行資料
x = data['x']
y = data['y']
print x
print y

#非線性最小二乘法擬合
popt, pcov = curve_fit(func, x, y)
#獲取popt裡面是擬合係數
a = popt[0] 
b = popt[1]
yvals = func(x,a,b) #擬合y值
print u'係數a:', a
print u'係數b:', b

#繪圖
plot1 = plt.plot(x, y, 's',label='original values')
plot2 = plt.plot(x, yvals, 'r',label='polyfit values')
plt.xlabel('x')
plt.ylabel('y')
plt.legend(loc=4) #指定legend的位置右下角
plt.title('curve_fit')
plt.savefig('test3.png')
plt.show()

4.三個引數擬合

        最後介紹官方給出的例項，講述傳遞三個引數，通常為 a*e(b/x)+c形式。

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

def func(x, a, b, c):
    return a * np.exp(-b * x) + c

# define the data to be fit with some noise
xdata = np.linspace(0, 4, 50)
y = func(xdata, 2.5, 1.3, 0.5)
y_noise = 0.2 * np.random.normal(size=xdata.size)
ydata = y + y_noise
plt.plot(xdata, ydata, 'b-', label='data')

# Fit for the parameters a, b, c of the function `func`
popt, pcov = curve_fit(func, xdata, ydata)
plt.plot(xdata, func(xdata, *popt), 'r-', label='fit')

# Constrain the optimization to the region of ``0 < a < 3``, ``0 < b < 2``
# and ``0 < c < 1``:
popt, pcov = curve_fit(func, xdata, ydata, bounds=(0, [3., 2., 1.]))
plt.plot(xdata, func(xdata, *popt), 'g--', label='fit-with-bounds')

plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.show()

三. 冪律分佈擬合及疑問

        下面是我冪率分佈的實驗，因為涉及到保密，所以只提出幾個問題。
        圖1是多項式的擬合結果，基本符合圖形趨勢。
        圖2是冪指數擬合結果，冪指數為-1.18也符合人類的基本活動規律。