效果
數(shù)據(jù)
來自動物學實驗-自選實驗
原理
邏輯斯蒂增長曲線是種群研究中一個非常經(jīng)典的增長曲線擒悬,它依據(jù)邏輯斯蒂方程(微分形式),描述了在食物、空間等資源有限時,由于種群數(shù)量增長醉冤,食物、空間資源逐漸進入“僧多粥少”的局面篙悯,從而導致種群增速放緩并最終進入一種動態(tài)平衡的狀態(tài)蚁阳。
負密度依賴調(diào)節(jié)下的種群增長可用邏輯斯蒂模型(logistic model)描述:
??:種群大小
??:內(nèi)稟增長率
??:環(huán)境承載力(carrying capacity)
- 1844年比利時數(shù)學家P. Verhulst提出;1921年美國生物學家R. Pearl和L. Reed再發(fā)現(xiàn)鸽照,用于預測人口增長
我們對這條式子進行變換螺捐,積分后得到:
進一步變換可以得到:
可以近似看成直線,內(nèi)稟增長率r等于斜率的相反數(shù)
python實現(xiàn)
拆解一下我們的任務:
1.準備數(shù)據(jù)集移宅,從數(shù)據(jù)里算得K值(環(huán)境容納量)
2.對數(shù)據(jù)進行變換归粉,使之結(jié)構(gòu)上可以看成的形式,也就是一條直線
3.把擬合的直線畫出來漏峰,再依據(jù)截距(intercept)和斜率(slope)反解出原方程中的參數(shù)r(內(nèi)稟增長率)和截距a
4.通過反解的參數(shù)把邏輯斯蒂增長曲線畫出來
導入numpy糠悼、matplotlib.pyplot、stats三個模塊
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
設置字體(如果你也打算用中文標簽浅乔,最好用全局字體倔喂,不然會后悔)
plt.rcParams['font.family'] = 'SimSun'
準備數(shù)據(jù)集铝条,量比較小,可以使用numpy庫下的array數(shù)組
有人就要問了:為什么不用pandas庫下面的Series和dataframe席噩?
好問題班缰,當時的我并不熟悉dataframe,現(xiàn)在的我懶得改( ??? )
想畫幾個就放幾個數(shù)據(jù)集進來悼枢,也可以單獨繪制某一個
name1 = '對照組'
t1 = np.array([1, 2, 3, 4, 5, 6, 7, 8])
N1 = np.array([0.5, 8, 17, 269, 690, 816, 1582, 1316])
K1 = 1322
name2 = 'ph = 6'
t2 = np.array([1, 2, 3, 4, 5, 6, 7, 8])
N2 = np.array([0.5, 8, 27, 194, 500, 670, 759, 654])
K2 = 694
name3 = 'ph = 8'
t3 = np.array([1, 2, 3, 4, 5, 6, 7, 8])
N3 = np.array([0.5, 12, 47, 424, 516, 468, 734, 910])
K3 = 846
最最最困難的地方埠忘,沒有沒有現(xiàn)成的函數(shù),需要自己定義一個:
我們希望這個函數(shù)能利用五個參數(shù)(t時間, N種群數(shù)量, K環(huán)境承載力, name組別, color顏色)馒索,實現(xiàn)后三個目標:
1.反解參數(shù)r
2.畫變換后數(shù)據(jù)的擬合直線
3.依據(jù)參數(shù)畫邏輯斯蒂曲線
stats.linregress可以反解我們需要的值莹妒,返回五個參數(shù),分別是斜率绰上、截距旨怠、相關系數(shù)、p值蜈块、斜率的標準誤差(此處用不到p值鉴腻,會發(fā)現(xiàn)p值全為零);
matplotlib.pyplot簡稱plt下有各種函數(shù)可供使用百揭,subplot為子圖分區(qū)爽哎、plot繪圖、scatter標注數(shù)據(jù)點信峻、xlabel和ylabel標注坐標軸倦青、title賦予圖片標題
def logistic_fit(t, N, K, name, color):
# 反解參數(shù)
slope, intercept, r_value, p_value, std_err = stats.linregress(t, np.log(abs(K - N) / N))
r = -slope
a = intercept
t_fit = np.linspace(min(t), max(t), 1000)
N_fit = K / (1 + np.exp(-r * t_fit + a))
print(f"{name}的內(nèi)稟增長率是{"%.3f" % r}")
print(f"{name}的截距是{"%.3f" % intercept}")
plt.subplot(1, 2, 1)
plt.plot(t, slope * t + intercept, color=color, label=name)
plt.scatter(t, np.log(abs(K - N) / N), label=name, color=color)
plt.xlabel('培養(yǎng)時間(天)')
plt.ylabel('變換后數(shù)據(jù):log((K - N) / N)')
plt.title('變換后數(shù)據(jù)的線性回歸')
plt.legend()
plt.subplot(1, 2, 2)
plt.plot(t_fit, N_fit, color=color, label=name)
plt.scatter(t, N, label=name, color=color)
plt.xlabel('培養(yǎng)時間(天)')
plt.ylabel('種群密度(單位:只/ml)')
plt.title('邏輯斯蒂回歸曲線')
plt.legend()
return t_fit, N_fit, slope, intercept
函數(shù)封裝好后,調(diào)用函數(shù)就好啦
t_fit1, N_fit1, slope1, intercept1 = logistic_fit(t1, N1, K1, name1, 'red')
t_fit2, N_fit2, slope2, intercept2 = logistic_fit(t2, N2, K2, name2, 'green')
t_fit3, N_fit3, slope3, intercept3 = logistic_fit(t3, N3, K3, name3, 'blue')
還有一件事盹舞!如果在函數(shù)中寫plt.show()就會導致所有圖片不在同一圖層产镐,而且是一張一張顯示,只需要把它拿開踢步,最后再進行統(tǒng)一的繪圖就可以了
plt.show()