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

2024-11-24 本文已影响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()