效果
效果图.png数据
来自动物学实验-自选实验
image.png原理
逻辑斯蒂增长曲线是种群研究中一个非常经典的增长曲线,它依据逻辑斯蒂方程(微分形式),描述了在食物、空间等资源有限时,由于种群数量增长,食物、空间资源逐渐进入“僧多粥少”的局面,从而导致种群增速放缓并最终进入一种动态平衡的状态。
负密度依赖调节下的种群增长可用逻辑斯蒂模型(logistic model)描述:
𝑁:种群大小
𝑟:内禀增长率
𝐾:环境承载力(carrying capacity)
- 1844年比利时数学家P. Verhulst提出;1921年美国生物学家R. Pearl和L. Reed再发现,用于预测人口增长
我们对这条式子进行变换,积分后得到:
进一步变换可以得到:
可以近似看成直线,内禀增长率r等于斜率的相反数
python实现
拆解一下我们的任务:
1.准备数据集,从数据里算得K值(环境容纳量)
2.对数据进行变换,使之结构上可以看成的形式,也就是一条直线
3.把拟合的直线画出来,再依据截距(intercept)和斜率(slope)反解出原方程中的参数r(内禀增长率)和截距a
4.通过反解的参数把逻辑斯蒂增长曲线画出来
导入numpy、matplotlib.pyplot、stats三个模块
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
设置字体(如果你也打算用中文标签,最好用全局字体,不然会后悔)
plt.rcParams['font.family'] = 'SimSun'
准备数据集,量比较小,可以使用numpy库下的array数组
有人就要问了:为什么不用pandas库下面的Series和dataframe?
好问题,当时的我并不熟悉dataframe,现在的我懒得改( ◠‿◠ )
想画几个就放几个数据集进来,也可以单独绘制某一个
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
最最最困难的地方,没有没有现成的函数,需要自己定义一个:
我们希望这个函数能利用五个参数(t时间, N种群数量, K环境承载力, name组别, color颜色),实现后三个目标:
1.反解参数r
2.画变换后数据的拟合直线
3.依据参数画逻辑斯蒂曲线
stats.linregress可以反解我们需要的值,返回五个参数,分别是斜率、截距、相关系数、p值、斜率的标准误差(此处用不到p值,会发现p值全为零);
matplotlib.pyplot简称plt下有各种函数可供使用,subplot为子图分区、plot绘图、scatter标注数据点、xlabel和ylabel标注坐标轴、title赋予图片标题
def logistic_fit(t, N, K, name, color):
# 反解参数
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}的内禀增长率是{"%.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('培养时间(天)')
plt.ylabel('变换后数据:log((K - N) / N)')
plt.title('变换后数据的线性回归')
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('培养时间(天)')
plt.ylabel('种群密度(单位:只/ml)')
plt.title('逻辑斯蒂回归曲线')
plt.legend()
return t_fit, N_fit, slope, intercept
函数封装好后,调用函数就好啦
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')
还有一件事!如果在函数中写plt.show()就会导致所有图片不在同一图层,而且是一张一张显示,只需要把它拿开,最后再进行统一的绘图就可以了
plt.show()
网友评论