美文网首页大数据 爬虫Python AI Sql数据科学与R语言
python绘图-草履虫种群增长(逻辑斯蒂方程)

python绘图-草履虫种群增长(逻辑斯蒂方程)

作者: 古代怪兽哥莫拉超进化 | 来源:发表于2024-11-24 23:13 被阅读0次

效果

效果图.png

数据

来自动物学实验-自选实验
image.png

原理

逻辑斯蒂增长曲线是种群研究中一个非常经典的增长曲线,它依据逻辑斯蒂方程(微分形式),描述了在食物、空间等资源有限时,由于种群数量增长,食物、空间资源逐渐进入“僧多粥少”的局面,从而导致种群增速放缓并最终进入一种动态平衡的状态。

负密度依赖调节下的种群增长可用逻辑斯蒂模型(logistic model)描述:
\frac{dN}{dt} = rN(1-\frac{N}{K})
𝑁:种群大小
𝑟:内禀增长率
𝐾:环境承载力(carrying capacity)

  • 1844年比利时数学家P. Verhulst提出;1921年美国生物学家R. Pearl和L. Reed再发现,用于预测人口增长
    我们对这条式子进行变换,积分后得到:
    N = K/(1+e^{a - rt})
    进一步变换可以得到:
    a - rt = ln[(K-N)/N]
    可以近似看成直线,内禀增长率r等于斜率的相反数

python实现

拆解一下我们的任务:

1.准备数据集,从数据里算得K值(环境容纳量)
2.对数据进行变换,使之结构上可以看成y = a + bx的形式,也就是一条直线
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()

相关文章

网友评论

    本文标题:python绘图-草履虫种群增长(逻辑斯蒂方程)

    本文链接:https://www.haomeiwen.com/subject/sowzjjtx.html