r语言 plot_R语言SIR模型（Susceptible Infected Recovered Model）代码实例

最新推荐文章于 2024-05-16 08:30:00 发布

weixin_39672680

最新推荐文章于 2024-05-16 08:30:00 发布

阅读量338

点赞数

文章标签： r语言 plot

本文链接：https://blog.csdn.net/weixin_39672680/article/details/111610044

版权

本文介绍了如何使用R语言实现SIR模型，该模型用于模拟信息传播。首先，文章定义了SIR模型的三个状态：易感者（S）、感染者（I）和移除者（R）。接着，详细阐述了模型的应用和微分方程，并展示了如何生成不同类型的网络。然后，通过设置传染和恢复概率，模拟了网络中的扩散过程。最后，通过可视化展示了SIR模型的扩散路径和时间演变。

摘要由CSDN通过智能技术生成

原文链接：

拓端数据科技 / Welcome to tecdattecdat.cn

SIR模型定义

SIR模型是一种传播模型，是信息传播过程的抽象描述。
SIR模型是传染病模型中最经典的模型，其中S表示易感者，I表示感染者，R表示移除者。

S：Susceptible，易感者I：Infective，感染者R：Removal，移除者

SIR模型的应用

SIR模型应用于信息传播的研究。

传播过程大致如下：最初，所有的节点都处于易感染状态。然后，部分节点接触到信息后，变成感染状态，这些感染状态的节点试着去感染其他易感染状态的节点，或者进入恢复状态。感染一个节点即传递信息或者对某事的态度。恢复状态，即免疫，处于恢复状态的节点不再参与信息的传播。

SIR的微分方程

a为感染率、b恢复率

注意：

t为某个时刻，例如t=1，S(1)为第一天易感人群的人数。
无论t为什么时刻，总人数是不变的，即N(t)=S(t)+I(t)+R(t)。
人口总数总保持一个常数，即N(t)=K，不考虑人口的出生、死亡、迁移等因素。

这里介绍一个使用R模拟网络扩散的例子。

第一步，生成网络。

规则网

g =graph.tree(size, children =2); plot(g)

g =graph.star(size); plot(g)

g =graph.full(size); plot(g)

g =graph.ring(size); plot(g)

第二步，随机选取一个或n个随机种子。

# initiate the diffusers
seeds_num =1 diffusers =sample(V(g),seeds_num) ;
diffusers
## + 1/50 vertex:
## [1] 43
infected =list()
infected[[1]]=diffusers#

第三步，传染能力

在这个简单的例子中，每个节点的传染能力是0.5，即与其相连的节点以0.5的概率被其感染,每个节点的回复能力是0.5，即其以0.5的概率被其回复。在R中的实现是通过抛硬币的方式来实现的。

## [1] 0

显然，这很容易扩展到更一般的情况，比如节点的平均感染能力是0.128，那么可以这么写：节点的平均回复能力是0.1，那么可以这么写

p =0.128
coins =c(rep(1, p*1000), rep(0,(1-p)*1000))
sample(coins, 1, replace=TRUE, prob=rep(1/n, n))
## [1] 0
n =length(coins2)
sample(coins2, 1, replace=TRUE, prob=rep(1/n, n))
## [1] 0

当然最重要的一步是要能按照“时间”更新网络节点被感染的信息。

keep =unlist(lapply(nearest_neighbors[,2], toss))
new_infected =as.numeric(as.character(nearest_neighbors[,1][keep >=1]))
diffusers =unique(c(as.numeric(diffusers), new_infected))
return(diffusers)}
set.seed(1);

开启扩散过程！

先看看S曲线吧：

为了可视化这个扩散的过程，我们用红色来标记被感染者。

# generate a palette#
plot(g, layout =layout.old)
set.seed(1)#
library(animation)# start the plot
m =1

如同在Netlogo里一样，我们可以把网络扩散与增长曲线同时展示出来：

set.seed(1)
# start the plot
m =1
p_cum=numeric(0)
h_cum=numeric(0)
i_cum=numeric(0)
while( m<50 ) {# start the plot
layout(matrix(c(1, 2, 1, 3), 2,2, byrow =TRUE), widths=c(3,1), heights=c(1, 1))
V(g)$color = "white"
V(g)$color[V(g)%in%infected[[m ]] ] = "red"
V(g)$color[V(g)%in%health[[m ]]] = "green"
if(m<=length(infected))
plot(pp~time, type ="h", ylab ="PDF", xlab ="Time",xlim =c(0,i), ylim =c(0,1), frame.plot =FALSE)
m =m +1
}