• Posts tagged "R"

Blog Archives

用R编写markdown格式文档

R的极客理想系列文章,涵盖了R的思想,使用,工具,创新等的一系列要点,以我个人的学习和体验去诠释R的强大。

R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

要成为有理想的极客,我们不能停留在语法上,要掌握牢固的数学,概率,统计知识,同时还要有创新精神,把R语言发挥到各个领域。让我们一起动起来吧,开始R的极客理想。

关于作者:

  • 张丹,分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-markdown-rmarkdown/

前言

Markdown格式文档编写,已经从小众成为了一种主流的文档编写方式。R语言社区在可重复性文档编写领域中,有着重要的贡献。knitr和rmarkdown等包的支持,到后来在RStudio中直接集成的rmarkdown的功能,让R的广大用户有了非常舒服的markdown文档使用体验。

2015时,我就已经开始用markdown来写书稿,当时写了2篇关于在nodejs领域中使用markdown的文章,Marked高效的Markdown解析器用WebStorm编辑Markdown

时光证明了,小而美的好东西,真的会放大光芒。改变用户使用习惯,从小而美开始。

目录

  1. 创建一个rmarkdown文档
  2. 嵌入markdown文本
  3. 嵌入R代码
  4. 嵌入DT动态表格
  5. 嵌入echarts动态图
  6. 嵌入数学公式

1. 创建一个rmarkdown文档

写文档的常规操作是用word,程序员写文档的常规操作是用IDE的编程工具,用markdown格式来写文档,最后可生成word、pdf、html等文档。

我们今天讲讲怎么用R语言来写文档吧,那么R语言对markdown的支持包是rmarkdown包,现在新版的RStudio开发工具,已经集成了rmarkdown包,我们直接在RStudio中操作就行了。

首先我们先打开RStudio开发工具,选择新建一个R Markdown文件。

选择Document文档格式,设置Title, Author, Date等基本信息,然后可以选择文档输出格式HTML, PDF, Word的三种选择,不同的格式输出需求本机已经安装对应的其他软件。

打开新创建的文档,默认会是一个模板,以rmarkdown的格式进行展示,上面所输入的参数,在文档都的开头部分都有展示。然后保存一下文件。

我们可以点击Knit的下拉菜单,选择输出为HTML。

然后,运行Knit就会新打开一个窗口,看到HTML的网页输出。

如果想输出为word格式,重新选择 Knit到word就行了,运行后会弹出word的新生成的文件。

一通操作下来,功能切换很是顺畅,体验非常舒服。那么接下来就在rmarkdown中,加入各种元素,让整个文档丰富起来。

2. 嵌入markdown文本

在rmarkdown中嵌入markdown文本,是最基本的操作,在默认生成的文档中,就是一个带有markdown格式的文本。markdown格式的文本,包括10种主要的格式,分别是标题,字体,引用,分隔,图片,链接,列表,表格,代码。

2.1 标题
标题支持六级标题,用#做为开头,#的数量表示是几级标题。


# 这是一级标题
## 这是二级标题
### 这是三级标题
#### 这是四级标题
##### 这是五级标题
###### 这是六级标题


图中,左边中编辑界面,右边是HTML的输出界面。

2.2 字体
字体格式设置有2种方法,一种是就markdown默认支持的4种格式,加粗,斜体,斜体加粗,删除线。另一种方法就是使用html的标签<font>来进行设置。


**这是加粗的文字**
*这是倾斜的文字*`
***这是斜体加粗的文字***
~~这是加删除线的文字~~

<font face="宋体">宋体</font>
<font face="宋体" color="red">红色宋体</font>
<font face="宋体" color="red" size="5">5号红色宋体</font>

2.3 引用
使用>用于表示引用外部的一段话,在引用中可以用\来进行换行。


> 我要引用一段文字\
\
R语言作为统计学一门语言,一直在小众领域闪耀着光芒。直到大数据的爆发,R语言变成了一门炙手可热的数据分析的利器。随着越来越多的工程背景的人的加入,R语言的社区在迅速扩大成长。现在已不仅仅是统计领域,教育,银行,电商,互联网….都在使用R语言。

2.4 分隔

在markdown中,可以通过3个连续的-或*,来设置分隔符。


--开始

---

**开始

********************

结束

2.5 图片
插入图片时,需要使用![文本说明](链接)语法。


![R的极客理想](http://fens.me/images/slider/img4.png)


2.6 链接
设置超级链接,有2种写法,一种就是[显示文字](链接),另一种是<链接>。


[粉丝日志博客](fens.me)

<http://www.fens.me>

2.7 列表
列表支持3种语法,一种是无序列表,二是有序列表,三是嵌套列表。

  • 无序列表,用 – + * 任何一种都可以,跟内容之间都要有一个空格
  • 有序列表,用数字加点,跟内容之间都要有一个空格
  • 嵌套列表,上一级和下一级之间有三个空格,跟内容之间都要有一个空格

无序列表

+ 列表内容
+ 列表内容
+ 列表内容

有序列表

1. 列表内容
2. 列表内容
3. 列表内容

无序列表嵌套

+ 1 列表内容
   + 1.1 列表内容
   + 1.2 列表内容
+ 2 列表内容
+ 3 列表内容

2.8 表格
在嵌入表格时,需求包括表头行,分隔符行,表体行。表头行用|分隔,表体行也用|分隔,分隔符行需要用–|组合进行分隔。当希望单元格整列对齐时,可以分隔符行加:来控制对齐,文字默认居左,两边加表示文字居中,右边加:表示文字居右。


表头|表头|表头
---|:--:|---:
内容|内容|内容
内容|内容|内容

表头|表头|表头
--|--|--
内容|内容|内容
内容|内容|内容


2.9 代码
代码的嵌入可以分为行内代码嵌入代码和代码块。

  • 行内代码,需要使用”包住代码,代码不可换行。
  • 代码块,需要用“` “`来包住代码,代码可换行。

文字中内嵌了`a="abc"`代码

```
a="abc"
a1=2
a2=3
a3=a1+a2
```

Markdown的基本语法,就是这9部分很简单,写一次练习一下,就记住了。

3. 嵌入R代码

上面我们演示了在markdown中嵌入代码,那么对于R语言来说,如果我们希望嵌入R代码,需要代码块针对R代码做一些支持的时候,要怎么写呢。


```{r}
a<-Sys.Date()
a

head(iris)
plot(iris)
```

接下来,我们还可以在{r echo=FALSE}加入echo=FALSE,这样就不输出R的代码本身了。


```{r echo=FALSE}
library(xts)
xdf<-xts(rnorm(100),order.by = as.Date('2022-01-01')+1:100)
plot(xdf)
```

4. 嵌入DT动态表格

在R的代码中,我们可以调用R的各种包进行输出,可以使用DT包来做动态表格,只要把数据放到datatable()函数中,就可以在markdown中生成对应的动态表格,支持翻页、关键字搜索,排序等多种功能。


```{r}
library(DT)
datatable(iris)
```

5. 嵌入echarts动态图

在文档输出时,经常需要用图进行展示,我们也可以直接使用echart4r包来调用echarts进行多种图的可视化展示。


```{r}
library(echarts4r)
library(magrittr)
iris %>% 
  e_charts(x = Sepal.Length) %>%
  e_scatter(serie = Sepal.Width )
```

6. 嵌入数学公式

文档中,也避免不了经常会写一些数据公式 ,rmarkdown中也集成是嵌入数据公式的功能。我们可以直接使用$$公式$$的语法,来生成复杂的数学公式。在文档中加入数学公式,是使用mathjax格式来进行数学公式的编辑。


$$
\hat\beta_1=\frac{\sum\left(X_i-\overline{X}\right)Y_i}{\sum\left(X_i-\overline{X}\right)^2}=\frac{\sum X_iY_i-\overline{X}\sum Y_i}{\sum X_i^2+\sum \overline{X}^2 - 2\overline{X}\sum X_i}\\
=\frac{\sum X_iY_i-\overline{X}*N\overline{Y}}{\sum X_i^2+N\overline{X}^2 - 2\overline{X}*N\overline{X}}
=\frac{\sum^{N}_{i=1} X_iY_i-N\overline{XY}}{\sum^{N}_{i=1}Xi^2-N\overline{X}^2}
$$

本文我介绍了rmarkdown的功能,有非常丰富的支持,多样的表达形式,语法简单,而且是可重复性生成的,比起 word 写文档,有更多的编程的特点。写书,写论文,写标书,写文档,都是不错的选择。

转载请注明出处:
http://blog.fens.me/r-markdown-rmarkdown/

打赏作者

知识图谱:最短路径算法Bellman-Ford

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-bellman-ford/

前言

最短路径是图算法中的基本算法,被普遍应用于多种场景中,如企业生产经营的场景中,现金流的生产和交换的过程,会产生正向的现金流比如经营收入所得,也会产生负向的现金流比如借款、支出等。如果我们以企业作为节点,以现金流做为边时,由于出现了负向的现金流的边,这样构成负权重的网络图谱,就需要用Bellman-Ford算法来实现。

常用的最短路径算法包括:Dijkstra算法,A star算法,Bellman-Ford算法,SPFA算法(Bellman-Ford算法的改进版本),Floyd-Warshall算法,Johnson算法以及Bi-direction BFS算法。

目录

  1. Bellman-Ford算法介绍
  2. R语言调用Bellman-Ford算法
  3. 自己手写Bellman-Ford算法
  4. 算一个复杂点的网络关系图

1. Bellman-Ford算法介绍

Bellman Ford算法最初于1955年被Alfonso Shimbel提出,但最终基于1958和1956年发表论文的 Richard Bellman和Lester Ford, Jr.二人命名。Edward F. Moore于1957年也提出了该算法,因此该算法有时也称作是Bellman Ford Moore算法。

Bellman-Ford算法与Dijkstra算法一样,解决的是单源最短路径问题。两者不同之处在于,Dijkstra只适用于无负权边的图,而Bellman-Ford算法无此限制:无论是否有向,无论是否有负权边,只要图中没有负权环,则该算法可以正确地给出起点到其余各点的最短路径,否则报告负权环的存在。

Bellman-Ford 算法采用动态规划进行设计,实现的时间复杂度为 O(V*E),其中 V 为顶点数量,E 为边的数量。每次松弛操作实际上是对相邻节点的访问(相当于广度优先搜索),第n次松弛操作保证了所有深度为n的路径最短。由于图的最短路径最长不会经过超过|V| – 1条边,所以可知贝尔曼-福特算法所得为最短路径,也可只时间复杂度为O(VE)。

Bellman-Ford 算法描述:

  • 创建源顶点 v 到图中所有顶点的距离的集合 distSet,为图中的所有顶点指定一个距离值,初始均为 Infinite,源顶点距离为 0;
  • 计算最短路径,执行 V – 1 次遍历;
  • 对于图中的每条边:如果起点 u 的距离 d 加上边的权值 w 小于终点 v 的距离 d,则更新终点 v 的距离值 d;
  • 检测图中是否有负权边形成了环,遍历图中的所有边,计算 u 至 v 的距离,如果对于 v 存在更小的距离,则说明存在环;

下面,我们来举例说明一下算法,创建一个带负权重的有向图,来描述对最短路径

开始点为A,数组U为开始列,数组S为上层节点列。通过6次遍历,找到从A到G的最短路径。

通过4次遍历,其实就能获得从A到D的最短路径。

2. R语言调用Bellman-Ford算法

接下来,我们可以继续使用R语言的igraph包来进行实现。新建数据集,如上面例子所示。

# 加载igraph库
> library(igraph)

# 新建图数据
> df<-data.frame(
+   from=c("A","A","B","B","B","D","D","E"),
+   to=c("B","C","C","D","E","B","C","D"),
+   weight=c(-1,4,3,2,2,1,5,-3)
+ )
 
# 用igraph画图 
> a<-graph_from_data_frame(df, directed = TRUE)
> E(a)$arrow.size<-1
> E(a)$arrow.width<-0.2
> E(a)$label<-df$weight

# layout需要设置一下,不然图的展示会很不好看
> plot(a,layout=layout.gem)

layout需要设置一下,不然图的展示会节点叠加在一起,就看不出来图的样子了,我这里选择layout.gem类型。

然后,使用shortest.paths()函数,进行bellman-ford算法计算最短路径的时候,就报错了。错误提示,图形成负权重的环路。

# 使用bellman-ford算法
> shortest.paths(a,"A", algorithm = "bellman-ford")
Error in shortest.paths(a, "A", algorithm = "bellman-ford") : 
  At core/paths/bellman_ford.c:157 : cannot run Bellman-Ford algorithm, Negative loop detected while calculating shortest paths

尝试了多种数据集,只要是有负值shortest.paths()函数都会报错,无论是否有负权重的环路,我也尝试找了几个的其他的R包,竟然没有发现Bellman-Ford算法的实现,因此只好自己手动写一个了。

3. 自已手写Bellman-Ford算法

于是,我就自己实现bellman-ford算法函数。写了2个小时,好像不练习了,可能还是会有bug吧。

在实现过程中,我们没有做“环”的判断,可能会导致计算逻辑出现问题。我假设了节点的寻路的过程是排序好的,不会出现对于一节点的反复寻路,当然这个假设条件,在某些复杂的场景中是不靠谱的,所以大家在使用的我的代码的时候,当一个基本版参考和学习用,不要直接用于自己的生产环境。

# bellman-ford 最短距离算法
# @param df数据集
# @param start开始节点
bellman_ford<-function(df,start="A"){
  
  # 初始化结果数据集
  rs<-data.frame(from=start,to=start,weight=0)
  
  # 初始化结束节点
  finishnode<-c()
  
  # 初始化开始节点
  nodes<-rs
  
  # 当开始节点数量<结束节点数量时,继续循环计算
  while(length(finishnode)<length(rs$to)){
    
    # 过滤已计算过的结束节点
    if(length(finishnode)>0){
      nodes<-rs[-which(rs$to %in% finishnode),]  
    }
    
    # 对原始数据进行循环,遍历每一行
    for(i in 1:nrow(df)){
      row<-df[i,];row
      
      # 对结束节点进行循环,遍历每一个结束节点
      for(j in nodes$to){
        if(row$from==j){
          # 把之前结束节点的权重叠加本次计算中
          row$weight<-row$weight+rs$weight[which(rs$to==j)]
          # 把结束节点加入到结果数据集
          rs<-rbind(rs,row)
        } 
      }
    }
    
    # 保留最短的唯一结束节点
    rs2<-rs[order(rs$weight),]
    n<-which(duplicated(rs2$to))
    if(length(n)<0){
      rs2<-rs2[-n,]
      rs<-rs2[order(as.numeric(row.names(rs2))),]
    }
    
    # 更新已经完成的结束节点
    finishnode<-c(finishnode,nodes$to)
  }
  return(rs)
}

执行计算,计算A点其他所有节点的距离。

> rs<-bellman_ford(df,"A")
> rs
  from to weight
1    A  A      0
2    A  B     -1
3    B  C      2
5    B  E      1
8    E  D     -2

结果数据集中,to 表示从开始点A到节点的距离,每一行from 到 to 表示,最短路径的顺序。

# 计算2点间的最短路径
# @param rs 结果数据集
# @param start 开始节点
# @param end 结束节点
get_shortest_path<-function(rs,start,end){
  
  # 初始化节点
  q<-c(start)
  idx<-which(rs$from==start)
  
  # 判断有几个分支
  if(length(idx)>0){
    for(i in idx){
      # 递归计算
      q2<-get_shortest_path(rs,rs$to[i],end)
      if(tail(q2,1)==end){
        q<-c(q,q2)
      }
    }
  }
  
  # 查看队列状态
  # print(q)
  
  # 当找到结束节点,程序停止
  if(tail(q,1)==end){
    return(q)
  }
  return(q)
}

# 计算从开始节点到结束节点的路径
> path<-get_shortest_path(rs[-1,],"A","D")
> path
[1] "A" "B" "E" "D"

最后,我们把从A点到D的最短路径涂上颜色,来出来。

# 把路径合并成data.frame
> path2<-data.frame(from=path[-length(path)],to=path[-1])
# 找到最短路径的索引
> idx<-which(paste0(df$from,"-",df$to) %in% paste0(path2$from,"-",path2$to))

# 涂色
> df$color<-"gray"
> df$color[idx]<-"red"
> E(a)$color<-df$color
> plot(a,layout=layout.gem)

红色路径就表示从A到D点的最短路径,看上去结果也是对的。在实现过程中,我们没有做“环”的判断,可能会导致计算逻辑出现问题。自己用写了一段算法,用R语言写算法,感觉比起C或者JAVA还要是简单很多的,因为大量基础函数提供了很多的功能和数据结构支持,不用像底层语言什么结构都需要自己从头写起。

4. 算一个复杂点的网路关系图

生成新数据集,共11个节点。

> df2<-data.frame(
+   from=c("A","B","C","D","C","F","B","B","F","G","H","F","E"),
+   to=  c("B","C","D","C","E","E","F","G","G","H","I","I","J"),
+   weight=c(5,20,10,-5,75,25,30,60,5,-50,-10,50,100)
+ )

用igraph进行画图展示。

> a2<-graph_from_data_frame(df2, directed = TRUE)
> E(a2)$arrow.size<-1
> E(a2)$arrow.width<-0.2
> E(a2)$label<-df2$weight
> plot(a2,layout=layout.gem)

计算A开始到所有节点的最短路径。

> rs<-bellman_ford(df2,"A")
> rs
   from to weight
1     A  A      0
2     A  B      5
3     C  D     35
6     F  E     60
7     B  F     35
9     F  G     40
10    G  H    -10
11    H  I    -20
13    E  J    160
21    B  C     25

计算A到D的最短路径路线,并画科展示。

# 计算 A 到 D的最短路径
> path<-get_shortest_path(rs[-1,],"A","D")
> path
[1] "A" "B" "C" "D"

# 画图
> path2<-data.frame(from=path[-length(path)],to=path[-1])
> idx<-which(paste0(df2$from,"-",df2$to) %in% paste0(path2$from,"-",path2$to))
> df2$color<-"gray"
> df2$color[idx]<-"red"
> E(a2)$arrow.size<-1
> E(a2)$arrow.width<-0.2
> E(a2)$color<-df2$color
> plot(a2,layout=layout.gem)

计算A到I的最短路径路线,并画科展示。

# 计算 A 到 I的最短路径
> path<-get_shortest_path(rs[-1,],"A","I")
> path
[1] "A" "B" "F" "G" "H" "I"

# 画图
> path2<-data.frame(from=path[-length(path)],to=path[-1])
> idx<-which(paste0(df2$from,"-",df2$to) %in% paste0(path2$from,"-",path2$to))
> df2$color<-"gray"
> df2$color[idx]<-"red"
> E(a2)$arrow.size<-1
> E(a2)$arrow.width<-0.2
> E(a2)$color<-df2$color
> plot(a2,layout=layout.gem)

算法和数据结构还是很重要的,图论算法原来仅在计算领域有应用,随着知识图谱的应用普及,相信未来会有更广阔的应用前景,把原来藏在背后的知识,扩展到前台形成新的模型思路。

努力吧,曾经的少年。

本文算法代码已上传到github: https://github.com/bsspirit/graph/blob/main/Bellman-Ford.r

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-bellman-ford/
打赏作者

知识图谱:最短路径算法Dijkstra

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-dijkstra/

前言

最短路径是图算法中的基本算法,被普遍应用于多种场景中,比如在实际生活中的路径规划、地图导航,建立道路交通网、物流运输网络、计算机网络等,这时就可以使用最短路径算法。

常用的最短路径算法包括:Dijkstra算法,A star算法,Bellman-Ford算法,SPFA算法(Bellman-Ford算法的改进版本),Floyd-Warshall算法,Johnson算法以及Bi-direction BFS算法。

目录

  1. 最短路径算法介绍
  2. Dijkstra 算法原理
  3. 用R语言调用 Dijkstra 算法

1. 最短路径算法介绍

最短路径问题是图论研究中一个经典算法问题,旨在寻找图中两节点或单个节点到其他节点之间的最短路径。根据问题的不同,算法的具体形式包括:

  • 确定起点的最短路径问题,即给定起始节点,求该节点到其他剩余节点的最短路径,适合使用Dijkstra算法
  • 确定终点的最短路径问题,即给定终点,求其他节点到该终点的最短路径。在无向图中,该问题与确定起点的问题等价;在有向图中,问题等价于把所有路径方向反转的确定起点的问题;
  • 确定起点终点的最短路径问题,即给定起点和终点,求两节点之间的最短路径;
  • 全局最短路径问题,即求图中所有节点之间的最短路径,适合使用Floyd-Warshall算法。

多种最短路径算法

最短路径算法 描 述
迪杰斯特拉算法
Dijkstra
寻找某个特定顶点到其它所有顶点的最短路径,该算法要求所有路径的权值为非负数。
贝尔曼福特算法
Bellman-Ford
寻找某个特定顶点到其它所有顶点的最短路径,该算法允许路径的权值为负数。
弗洛伊德算法
(Floyd-Warshall)
寻找各个顶点之间的最短路径,允许非环路的路径权值为负数,该算法不仅适用于稀疏图,在稠密图(路径数量多的图)中寻找最短路径的效率也很高。
约翰逊算法
(Johnson)
寻找各个顶点之间的最短路径,允许非环路的路径权值为负数,该算法更适用于稀疏图(路径数量少的图)。

2. Dijkstra 算法原理

Dijkstra算法,于1956年由荷兰计算机科学家艾兹赫尔.戴克斯特拉提出,用于解决赋权有向图的单源最短路径问题。所谓单源最短路径问题是指确定起点,寻找该节点到图中任意节点的最短路径,算法可用于寻找两个城市中的最短路径或是解决著名的旅行商问题。

基本思想

  1. 通过Dijkstra计算图中的最短路径时,需要指定一个起点A。
  2. 简历两个数组S和U,S用于记录已求出最短路径的顶点(以及相应的最短路径长度),U用于记录还未求出最短路径的顶点(以及该顶点到起点A的距离)。
  3. 初始时,数组S中只有起点A;数组U中是除起点A之外的顶点,并且数组U中记录各顶点到起点A的距离。如果顶点与起点A不相邻,距离为无穷大INF。
  4. 从数组U中找出路径最短的点X,并将其加入到数组S中;同时,从数组U中移除顶点X。接着,更新数组U中的各顶点到起点A的距离。
  5. 重复第4步操作,直到遍历完所有顶点。

下面,我们来举例说明一下算法,创建一个带权重的有向图,来描述对最短路径。

开始点为A,数组U为开始列,数组S为上层节点列。通过6次遍历,找到从A到G的最短路径。

3. 用R语言调用 Dijkstra 算法

我们用R语言中,igraph包已经实现了Dijkstra算法,我们可以直接调用igraph包的shortest.pahs()函数进行最短路径的计算。首先,安装和加载 igraph 程序包

> install.pacakges("igraph")
> library(igraph)

生一个有向有权重的网络图,如上例所示。

> df<-data.frame(
+   from=c("A","A","A","B","D","B","C","D","C","F","E","F"),
+   to=c("B","D","C","C","C","E","E","F","F","E","G","G"),
+   weight=c(4,6,6,1,2,7,6,5,4,1,6,8)
+ )
> df
   from to weight
1     A  B      4
2     A  D      6
3     A  C      6
4     B  C      1
5     D  C      2
6     B  E      7
7     C  E      6
8     D  F      5
9     C  F      4
10    F  E      1
11    E  G      6
12    F  G      8

把数据生成igraph图对象,再以网络图进行展示。

# 生成igraph对象
> a<-graph_from_data_frame(df, directed = TRUE, vertices = NULL)

# 画图
> E(a)$arrow.size<-1
> E(a)$arrow.width<-0.2
> E(a)$label<-df$weight
> plot(a)

上面我们把节点涂色了,下面继续把边也涂色。首先,找到最短路径上的边,把最短路径边涂红色,非最短路径的边涂灰色。

# 通过最短路径的节点,建立边的关系
> s_vert<-data.frame(from=s_node[-length(s_node)],to=s_node[-1])

# 找到网络图中所有的边
> a_vert<-as_data_frame(a,what="edges")

# 找到最短路径的边的索引
> idx<-which(paste0(a_vert$from,"-",a_vert$to) %in% paste0(s_vert$from,"-",s_vert$to))
> idx
[1]  1  4  9 10 11

# 涂色:红色为最短路径的边,灰色为非最短路径的边
> a_vert$color<-"gray"
> a_vert$color[idx]<-"red"

# 画图
> E(a)$color<-a_vert$color
> plot(a)

用Dijkstra算法,计算从A开始,到B,C,D,E,F,G的所有节点的最短距离。

> shortest.paths(a,"A", algorithm = "dijkstra")
  A B D C F  E  G
A 0 4 6 5 9 10 16

输出结果,从A到B最短距离为4,从A到D最短距离为6,从A到C最短距离为5,从A到F最短距离为9,从A到E最短距离为10,从A到G最短距离为16。

当然,我们除了想知道最短距离,还想知道最短距离的路径是什么样的。接下来,我们先计算最短距离路径上的点,并发点进行涂色。

# 计算最短距离经过的点
> s_node<-all_shortest_paths(a,"A","G")$res[[1]]$name
> snode
[1] "A" "B" "C" "F" "E" "G"

# 把节点进行涂色
> cols<-data.frame(label=V(a)$name,color=0)
> cols$color[which(cols$label %in% s_node)]<-1
> cols
  label color
1     A     1
2     B     1
3     D     0
4     C     1
5     F     1
6     E     1
7     G     1

# 画图
> V(a)$color<-cols$color
> plot(a)

其中,黄色点为最短距离穿过的节点,白色节点为不在最短距离路径上的节点。

上面我们把点进行了涂色,下面我们继续边最少路径的边进行涂色。首先,根据最短路径的节点,创建边的关系。然后对边进行涂色,在最短路径的边涂为红色,不在最短路径的边涂为灰色。

# 建立最短路径节点的关系
> s_vert<-data.frame(from=s_node[-length(s_node)],to=s_node[-1])

# 找到网络图的完整路径
> a_vert<-as_data_frame(a,what="edges")

# 找到最短路径的边的索引
> idx<-which(paste0(a_vert$from,"-",a_vert$to) %in% paste0(s_vert$from,"-",s_vert$to))
> idx
[1]  1  4  9 10 11

# 涂色
> a_vert$color<-"gray"
> a_vert$color[idx]<-"red"

# 画图
> E(a)$color<-a_vert$color
> plot(a)

这样就在完整的网络图中,可以从开始节点A,按红色路径走到G,生成了最短路径所过过的节点和边了。

扩展一下,我们用Dijkstra算法,分别计算从不同点开始,到所有不同的点的最短距离,通过矩阵来表达,任意2个节点的最短距离。

> shortest.paths(a, algorithm = "dijkstra")
   A  B  D  C F  E  G
A  0  4  6  5 9 10 16
B  4  0  3  1 5  6 12
D  6  3  0  2 5  6 12
C  5  1  2  0 4  5 11
F  9  5  5  4 0  1  7
E 10  6  6  5 1  0  6
G 16 12 12 11 7  6  0

通过上面的步骤,我们就实现了在R语言中使用Dijkstra算法,计算的最短路径算法。本来还想着,自己再来实现一遍,发现网上有很多的开源的代码,还让自己休息一下吧。下篇文章将介绍,知识图谱:最短路径算法Bellman-Ford算法,应该是需求自己手写了。

本文算法代码已上传到github: https://github.com/bsspirit/graph/blob/main/dijkstra.r

转载请注明出处:
http://blog.fens.me/r-graph-shortest-path-dijkstra/
打赏作者

知识图谱:社区发现算法leiden

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-graph-leidan/

前言

基于知识图谱的技术,已经越来越普遍的应用于各种业务场景中了,比如互联网领域的搜索推荐、智能客服;金融领域的智能风控、反欺诈识别;水利领域的最近的数字孪生建设规划也加入了知识平台的建设。

图算法方向,有很多的算法模型,本文将介绍Leiden算法,用于社区检测,帮助我们了解大型复杂网络的结构。新的技术趋势已经打开,抓住时代的机遇,用技术提升社会生产力。

目录

  1. Leiden算法介绍
  2. 用R语言加载leidenAlg包

1. Leiden算法介绍

社区检测通常用于了解大型复杂网络的结构。Leiden 算法是一种在大型关系网络上的社区检测算法,leiden算法计算社区之间的节点和边的关系,保证社区内部连接的有效性,并实现社区的分隔。当迭代应用 Leiden 算法时,它收敛到一个分区,在该分区中所有社区的所有子集都被局部最优分配,从而产生保证连接的社区。

leidan算法的论文:https://www.nature.com/articles/s41598-019-41695-z

Leiden 算法是为了改进 Louvain 算法的缺陷,Louvain 算法可能会发现任意连接不良的社区。 因此,Leiden保证了社区之间的良好连接, 而且算法速度更快,可扩展性更好,并且可以在数百万个节点的图上运行(只要它们可以放入内存)。

Leiden 算法包括三个阶段:

  1. 节点局部移动,快速找到分区
  2. 分区细化
  3. 基于细化分区的网络聚合,利用非细化分区为聚合网络创建初始分区。 迭代步骤直到收敛。

2. leidenAlg包算法实现

R语言中,已经有了Leiden算法的实现,Leiden算法是在2019年刚推出,R社区的开发效率真是很高!我们可以加载leidanAlg包来使用Leiden算法,结合igraph进行社区发现可视化展示。

项目地址:https://github.com/kharchenkolab/leidenAlg

首先,安装和加载leidanAlg包。安装过程很简单,如果igraph之前装过,可能会需要升级igraph到高版本。

# 安装leidenAlg 和 igraph包
> install.packages("leidenAlg")
> install.packages("igraph")

# 加载程序包
> library(leidenAlg)
> library(igraph)

# 设置工作目录
> setwd("C:/work/R/graph")

2.1 社区发现算法

首先,我先使用leidenAlg提供的一个数据集exampleGraph,进行算法的测试。

# 加载数据集exampleGraph
> data(exampleGraph)

# 查看数据集数据
> exampleGraph
IGRAPH 05dba3a UNW- 100 242 -- 
+ attr: name (v/c), weight (e/n), type (e/n)
+ edges from 05dba3a (vertex names):
 [1] MantonBM1_HiSeq_1-GGAACTTCACTGTCGG-1--MantonBM2_HiSeq_1-CTGATAGAGCGTTCCG-1
 [2] MantonBM1_HiSeq_1-ACGCCAGAGTACCGGA-1--MantonBM2_HiSeq_1-GGAAAGCCAGACGCCT-1
 [3] MantonBM1_HiSeq_1-ATCACGAAGGAGTCTG-1--MantonBM2_HiSeq_1-TTTGCGCGTAGCGCTC-1
 [4] MantonBM2_HiSeq_1-TTTGCGCGTAGCGCTC-1--MantonBM1_HiSeq_1-TTGGCAATCTTGAGAC-1
 [5] MantonBM2_HiSeq_1-TTTGCGCGTAGCGCTC-1--MantonBM1_HiSeq_1-ACGGGTCTCCACGACG-1
 [6] MantonBM2_HiSeq_1-TTTGCGCGTAGCGCTC-1--MantonBM1_HiSeq_1-GTCATTTGTCAGATAA-1
 [7] MantonBM2_HiSeq_1-TTTGCGCGTAGCGCTC-1--MantonBM1_HiSeq_1-CAGGTGCTCCACGTGG-1
 [8] MantonBM2_HiSeq_1-TTTGCGCGTAGCGCTC-1--MantonBM1_HiSeq_1-TTGGCAAGTCAGAATA-1
+ ... omitted several edges

可视化展示数据集的网路关系

> e<-exampleGraph
# 不现实节点名
> V(e)$label<-""

# 画图
> plot(e)

从上图中,可以看到整个网络是由多个小的社区组成的,肉眼是可以明显进行划分的。接下来,通过leiden算法,进行社区自动化的分隔。

# leidan算法计算
> lc<-leiden.community(e)

# 把不同的社区涂上不同的颜色
> library(RColorBrewer)
> node.cols <- brewer.pal(max(c(3, lc$membership)),"Pastel1")[lc$membership]
Warning message:
In brewer.pal(max(c(3, lc$membership)), "Pastel1") :
  n too large, allowed maximum for palette Pastel1 is 9
Returning the palette you asked for with that many colors

# 画图
> plot(e, vertex.color = node.cols, vertex.label="")

涂色后,可以直接进行观察,每个社区都刚好形成独立的颜色。leiden算法,可以自动地进行社区的识别。

2.2 递归聚类

递归 leiden算法, 使用 rleiden.community()函数 构建一个 n 步递归聚类,并返回树型图,用来展示社区各个对象之间的聚类关系。

建立深度为2的社区树型结构。

> r2 = rleiden.community(e, n.cores=1,max.depth=2)
> rs2<-as.dendrogram.fakeCommunities(r2)
> plot(rs2)

2.3 权重设置

构建一个新的网络关系图谱,我们用于边的权重,对社区划分的计算的影响。

首先,建立10节点的图谱,并设置边的连接和方向。


# 建立节点
> g <- make_star(10) 
> g<-add_edges(g,c(2,3,3,4, 4,5,2,4,2,5))
> g<-add_edges(g,c(6,7,6,8,6,9))
> E(g)$arrow.size<-0.2
> E(g)$arrow.width<-0.2

# 查看图谱结构
> g
IGRAPH d3bc7a7 D--- 10 17 -- In-star
+ attr: name (g/c), mode (g/c), center (g/n), arrow.size (e/n), arrow.width (e/n)
+ edges from d3bc7a7:
 [1]  2->1  3->1  4->1  5->1  6->1  7->1  8->1  9->1 10->1  2->3  3->4  4->5  2->4  2->5  6->7  6->8  6->9

# 画图
> plot(g)

给边设置权重,分别把1,2,3以随机数的方式,设置到不同的边上。

 # 取随机权重,设置边的大小,设置边的颜色
> E(g)$width <- sample(1:3,ecount(g),replace = TRUE)
> E(g)$color<-E(g)$width
> E(g)$width 
 [1] 1 2 1 1 2 1 3 1 3 3 2 3 3 2 3 1 2

# 画图
> plot(g)

上图中,按照权重的关系,我们把边凃上了颜色,可以查看图中节点和边的值。


# 查看图中节点和边的值
> as_data_frame(g)
   from to arrow.size arrow.width width color
1     2  1        0.2         0.2     1     1
2     3  1        0.2         0.2     2     2
3     4  1        0.2         0.2     1     1
4     5  1        0.2         0.2     1     1
5     6  1        0.2         0.2     2     2
6     7  1        0.2         0.2     1     1
7     8  1        0.2         0.2     3     3
8     9  1        0.2         0.2     1     1
9    10  1        0.2         0.2     3     3
10    2  3        0.2         0.2     3     3
11    3  4        0.2         0.2     2     2
12    4  5        0.2         0.2     3     3
13    2  4        0.2         0.2     3     3
14    2  5        0.2         0.2     2     2
15    6  7        0.2         0.2     3     3
16    6  8        0.2         0.2     1     1
17    6  9        0.2         0.2     2     2

最后,用leiden算法进行对节点的社区划分并涂色,就能看到最后,被社区分隔的结果。

# 根据权重进行社区划分计算,设置节点的颜色
> V(g)$color<-find_partition(g, E(g)$width)

# 画图
> plot(g)

权重的设置,可以改进leiden算法的社区发现的识别结果,自动地帮我们把强关联社区进行锁定,可以提升智能化的社区发现效率。

好的技术,用在合适的场景中,事半功倍。

本文算法代码已上传到github: https://github.com/bsspirit/graph/blob/main/leiden.r

转载请注明出处:
http://blog.fens.me/r-graph-leidan/

打赏作者

用R模拟本次北京新冠疫情

算法为王系列文章,涵盖了计算机算法,数据挖掘(机器学习)算法,统计算法,金融算法等的多种跨学科算法组合。在大数据时代的背景下,算法已经成为了金字塔顶的明星。一个好的算法可以创造一个伟大帝国,就像Google。

算法为王的时代正式到来….

关于作者:

  • 张丹, 分析师/程序员/Quant: R,Java,Nodejs
  • blog: http://blog.fens.me
  • email: bsspirit@gmail.com

转载请注明出处:
http://blog.fens.me/r-coivd19-beijing

前言

新冠疫情已经持续三年了,前段时间在北京和上海两个超大城市爆发,一度让工作和生活陷入恐慌。当时居家办公,每天核酸,都已经变成了生活的常态。

北京疫情从4月底爆发到,到5月全面严格管控,再到6月复工复产,虽然我们都亲身经历了,但并不知道为什么北京的政策是有效的。回顾一下,让数据说话,看看我们到底我们经历了什么。

由于本人非传染病专业,文中关于传染病的专业内容大都从网上汇集和再整理,当出现本文描述与教课书不符时,请以专业教材为准。

目录

  1. 获取北京新冠疫情数据
  2. R0基本传染数介绍
  3. 用SIR模拟北京疫情,结果失真
  4. 用SIR拟合算出R0,结果合理

1. 获取北京新冠疫情数据

在全球抗击新冠肺炎疫情的过程中,传染病动力学模型发挥了巨大的作用。传染病动模型,根据传染病的传播速度不同,空间范围各异,传播途径多样,动力学机理等各种因素,对传染病模型按照传染病的类型划分为 SI,SIR,SIRS,SEIR 模型。

首先,安装nCov2019,加载nCov2019包,再加载数据。关于nCov2019包的具体安装和使用,请参考文章:用R语言获取全球新冠疫情数据

# 设置工作路径
> setwd("C:/work/R/covid19")

# 如果没安装,请先安装nCov2019包
> #remotes::install_github("YuLab-SMU/nCov2019")

# 加载nCov2019包
> library(nCov2019)

# 下载最近数据
> x <- query()
Querying the latest data...
last update: 2022-07-05 
Querying the global data...
Gloabl total  555324897  cases; and  6362663  deaths
Gloabl total affect country or areas: 230
Gloabl total recovered cases: 165929
last update: 2022-07-05 
Querying the historical data...
Querying the vaccine data...
Total Candidates Programs : 51 
Querying the therapeutics data...
Total Candidates Programs : 84 
Query finish, each time you can launch query() to reflash the data

加载其他数据处理R包,并定义公共画图函数。Y轴做了log10的转化处理。

> library(ggplot2)
> library(magrittr)
> library(reshape2)
> 
> # 画图函数
> draw<-function(df){
+   dat<-melt(df,id.vars = c("date"))
+   g<-ggplot(data = dat, mapping = aes(x=date,y=value,colour=variable))
+   g<-g+geom_line() + geom_point()   
+   g<-g+scale_y_log10()
+   g<-g+ggtitle("北京疫情统计")+xlab("日期")+ylab("Log(人数)")
+   g
+ }

获取中国北京的新冠数据。

> hist<-x$historical
> beijing<-hist["China","beijing"]

# 查看数据集
> head(beijing)
    country province       date cases deaths recovered
6     China  beijing 2020-01-22    14      0         0
95    China  beijing 2020-01-23    22      0         0
184   China  beijing 2020-01-24    36      0         1
273   China  beijing 2020-01-25    41      0         2
362   China  beijing 2020-01-26    68      0         2
451   China  beijing 2020-01-27    80      1         2

计算每日新增数据,并画图展示北京整体疫情的情况。

> beijing$daily<-c(0,diff(beijing$cases))

# 画图
> draw(beijing[,c("date", "cases" ,"deaths", "recovered","daily")])

上图中cases为确诊人数,deaths为死亡人数,recovered为治愈人数(治愈后不会再得),daily为每日新增确诊人数。x轴为时间日期,y轴为人数(log化的人数)。

本次北京的疫情,从4月开始爆发,到6月基本控制住了。我们截取取从4月1日到7月4日数据,进行分析。

# 取2022-04-01开始到2022-07-04数据
> bj202204<-beijing[which(beijing$date>=as.Date("2022-04-01")),]
> head(bj202204)
      country province       date cases deaths recovered daily
71206   China  beijing 2022-04-01  1777      9         0     8
71295   China  beijing 2022-04-02  1777      9         0     0
71384   China  beijing 2022-04-03  1781      9         0     4
71473   China  beijing 2022-04-04  1791      9         0    10
71562   China  beijing 2022-04-05  1795      9         0     4
71651   China  beijing 2022-04-06  1800      9         0     5

# 计算这段时间的累计确诊数量
> bj202204$cum<-cumsum(bj202204$daily) 

# 画图
> draw(bj202204[,c("date", "cum","daily")])

上图中, daily为每日确诊数量,cum为累计确诊数量。x轴为时间日期,y轴为人数(log化的人数)。

2. R0基本传染数和Re有效传染数

在有了数据后,我们接下来要先了解一下R0基本传染数的概念,后面在用传染病动力学模型做拟合时,就可以用来对实际的疫情进行评价。

2.1 R0基本传染数

R0基本传染数,这是一个很重要的传染病学概念,用来衡量病原体的传染性。R0是指在一个完全易感人群中(没有任何预防手段介入并且所有人对此病原体没有免疫力的情况下),一个病例能传染的平均人数。

R0最常用的方法是使用累积发病率数据。

R0 = 传染期 * 每人每天的接触者数 * 每个接触者的感染概率
  • R0 < 1,意味着在没有防控的情况下, 一个感染者将传给少于一个人,所以传染病将会自行逐渐消失。
  • R0 > 1,每个感染者能传染多于一个人,那么就会有发展成流行病(epidemic)的可能性。但是一般不会永远持续,因为可能被感染的人口会慢慢减少。部分人口可以死该传染病,部分则可能治愈后产生免疫力。
  • R0 = 1,传染力维持,病院持续存在,会变成地方性流行病。

但R0不代表当前被感染的总人数。它也不是疾病严重程度的衡量指标。它只表示一个感染者平均能感染几个人,而不是告诉我们这些感染者疾病的严重程度。

对于一个特定的传染病,R0不一定是恒定的,它可能取决于人口密度和接触方式等因素。尽管如此,对于不同人群,一个特定传染病的R0是相对相似的。

部分传染病的R0值,数据来源 https://zh.m.wikipedia.org/zh-cn/基本传染数

疾病 传播途径 R0
麻疹 空气传播 12-18
白喉 唾液 6-7
天花 空气传播、 飞沫传播 5-7
脊髓灰质炎 粪口传播 5-7
风疹 空气传播、 飞沫传播 5-7
流行性腮腺炎 呼吸道飞沫 10-12
HIV/AIDS 性传播 4.5
百日咳 空气传播、 飞沫传播 5.5
SARS 空气传播、 飞沫传播、粪口传播 2-5
流行性感冒(1918年流感大流行) 空气传播、 飞沫传播 2-3
埃博拉出血热(西非埃博拉病毒疫症) 体液传播 1.5-2.5
COVID-19 飞沫传播、接触传播、粪口传播 1.4-3.8
3.8–8.9(截至2020年6月)
Delta变异株:5.1
Alpha变异株:4–5
Omicron变异株:7
中东呼吸综合征 呼吸道飞沫 0.5 (0.3 –0.8 )
普通感冒 呼吸道飞沫 2–3
水痘 空气传播疾病 10–12

有了上面的不同传染病的比较级,就能对本次新冠疫情有了一个数量上的感觉,Omicron变异株:7的数值,确实是非常高的传染性。

2.2 Re有效传染数

相对于R0,另一个指标Re可以帮我们分析,如何抑制疫情的爆发。Re为有效传染数,是指一个群体在任何特定的时间,不完全易感或已采取适当防控措施的情况下,单例病例能感染的平均人数。随着人口越来越多地获得免疫,无论是通过感染或接种疫苗后的个体免疫,还是随着人们的死亡,它都会发生变化。

Re = R0 *(1 – 有效控制率)*(易感人群比例)
  • 当Re>1时,每个受感染的个体将疾病传播给不止一个人,疾病可能在人群中传播。
  • 当Re<1时,可能会发生个体病例传播给另一个人发生新的感染,要遏制疾病流行,Re必须小于1,疫苗通过减少易感人群的比例,是使Re < 1的最佳手段。

3. 用SIR模型模拟北京疫情

我们首先用SIR模型进行北京疫情的模拟,具体操作,使用确定的SIR公式,通过设置参数带入模型模拟本次疫情的情况,观察SIR模型模拟的数据是否与实际发生的确诊数据,是否是匹配的。

我们直接使用EpiModel包,自带的SIR模型进行拟合。关于EpiModel包的详细使用介绍,请参考文章 EpiModel快速上手传染病模型

R程序实现代码

> library(EpiModel)
> param <- param.dcm(inf.prob = 0.5, act.rate = 0.5, rec.rate = 1/7)
> init <- init.dcm(s.num = 30*1000*1000, i.num = 8, r.num = 0)
> control <- control.dcm(type = "SIR", nsteps = 200, dt = 0.5)
> mod <- dcm(param, init, control)

设置模型初始参数:日感染率inf.prob=0.5,日传播率 act.rate=0.5,恢复率rec.rate=1/7。易感染者s.num=30*1000*1000,为北京人口数量,已感染者i.num=0,治愈r.num=0,观察周期nsteps=200。

> mod
EpiModel Simulation
=======================
Model class: dcm

Simulation Summary
-----------------------
Model type: SIR
No. runs: 1
No. time steps: 200
No. groups: 1

Model Parameters
-----------------------
inf.prob = 0.5
act.rate = 0.5
rec.rate = 0.1428571

Model Output
-----------------------
Variables: s.num i.num r.num si.flow ir.flow num

模型模拟结果

> plot(mod)

把实际的北京确诊数据,与模拟的确认数据画到一起。


# 生成日期列
> moddf<-as.data.frame(mod)
> moddf$date<-as.Date("2022-04-01")+0:(nrow(moddf)-1)

# 合并数据集
> d1<-data.frame(date=moddf$date,i=moddf$i.num,variable="sir")
> d2<-data.frame(date=bj202204$date,i=bj202204$cum,variable="real")
> mdf<-rbind(d1,d2)

# 画图
> g<-ggplot(data = mdf, mapping = aes(x=date,y=i,colour=variable))
> g<-g+geom_line() + geom_point()   
> g<-g+scale_y_log10()
> g<-g+ggtitle("北京疫情统计")+xlab("日期")+ylab("Log(人数)")
> g

红色为北京实际发生确诊的数据,蓝色线为SIR模型模拟的数据,从图中看是有很大的差异性的,所以直接用确定性的公式做疫情的模拟是很不靠谱的。

4. 基于 SIR 拟合北京疫情

使用SIR模型进行拟合,并设置运行参数


# 加载程序包
> library(lubridate)
> library(plyr)

# 设置初始化参数
> cases <- bj202204$cum     # 每日累计
> N <- 30*1000*1000         # 总人数
> startDay <- 6             # 开始累计时间

# 计算每日累计确诊
> infected <- cases[startDay:length(cases)]
> infected
 [1]   31   40   48   52   53   57   57   59   62   63   63   70   73   75   76   77   85  107  122  154
[21]  185  234  283  333  386  422  473  520  563  618  664  708  741  803  828  863  906  938  971 1010
[41] 1053 1105 1155 1206 1264 1316 1399 1440 1481 1518 1541 1561 1575 1583 1599 1613 1613 1632 1641 1641
[61] 1663 1663 1667 1669 1678 1717 1751 1780 1822 1847 1862 1875 1877 1878 1882 1886 1891 1895 1900 1905
[81] 1906 1909 1910 1911 1914 1914 1914 1915 1916 1920

# 以多少天为窗口计算R0
> window <- 12

定义SIR模型的拟合函数,用来进行后面的数据拟合。SIR模型的计算公式详解,请参考文章:用R语言解读传染病模型


# 定义SIR模型的拟合函数
> SIR<-function(time,vars,params){
+   with(as.list(c(vars,params)),{ #beta=rio,gamma=mu
+     num <- s.num + i.num + r.num
+     dS<- -beta * i.num * s.num/num
+     dI<- beta * i.num * s.num/num - gamma * i.num
+     dR<- gamma * i.num
+     return(list(c(dS,dI,dR)))
+   })
+ }

这里没有直接使用EpiModel包的SIR函数的原因是,EpiModel的参数命名规则不太一样,另外需要设置额外的参数,在这里用不上,比如inf.prob、act.rate等。当然如果一定用EpiModel的SIR函数,可以使用mod_SI_1g_cl(t, t0, parms)函数,自己做一下适配,并不复杂。


# 查看mod_SIR_1g_cl函数定义
> mod_SIR_1g_cl
function (t, t0, parms) 
{
    with(as.list(c(t0, parms)), {
        num <- s.num + i.num + r.num
        lambda <- inf.prob * act.rate * i.num/num
        if (!is.null(parms$inter.eff) && t >= inter.start) {
            lambda <- lambda * (1 - inter.eff)
        }
        si.flow <- lambda * s.num
        ir.flow <- rec.rate * i.num
        dS <- -si.flow
        dI <- si.flow - ir.flow
        dR <- ir.flow
        list(c(dS, dI, dR, si.flow, ir.flow), num = num)
    })
}

接下来,定义损失函数,由于后续采用的是最小二乘,所以都是适用与回归的损失函数。


> mse <- function(infected, fit) {
+   sum((infected - fit)^2) / length(infected)
+ }
> r2 <- function(infected, fit) {
+   1 - mse(infected, fit)/var(infected)
+ }

定义R0的拟合函数

> R0 <- function(days){
+ 
+   # 取窗口期内的数据拟合
+   WindowInfected <- infected[(days - window):(days - 1)]
+   # 持续天数
+   Days <- 1:(length(WindowInfected))
+   
+   # 设置SIR模型的默认值,s.num易感人群,i.num已感人群,r.num已治愈人群
+   init <- c(s.num = N - WindowInfected[1], i.num = WindowInfected[1], r.num = 0)
+ 
+   # 定义损失函数
+   LOSS <- function(parameters) {
+     names(parameters) <- c("beta", "gamma")
+     out <- ode(y = init, times = Days, func = SIR, parms = parameters)
+     fit <- out[ , 3]
+     -r2(WindowInfected,fit)
+   }
+ 
+   # 拟合计算
+   Opt = optim(c(0.5, 0.5), LOSS, method = "L-BFGS-B", lower = c(0, 0), upper = c(Inf, 1))
+   Opt_par = setNames(Opt$par, c("beta", "gamma"))
+   data.frame(days = days,
+              R0 = Opt_par["beta"] / Opt_par["gamma"],
+              beta = Opt_par["beta"],
+              gamma = Opt_par["gamma"],
+              value = Opt$value,
+              date = ymd("2022-04-01") + days(startDay + days - 1) - 1) # StartDay - 1 | days - 1
+ }

从13个日期2022-04-18开始,对每一个时间进行拟合,计算R0值。


> # 找到需要计算的日期
> days<-(window+1):(length(infected));days
 [1] 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46
[35] 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
[69] 81 82 83 84 85 86 87 88 89 90

> # 计算R0的计算
> RODF<-lapply(days,R0) %>% ldply
> head(RODF)
  days       R0      beta     gamma      value       date
1   13 1.163701 0.5813041 0.4995306 -0.5091804 2022-04-18
2   14 1.115407 0.5574453 0.4997684 -0.7966765 2022-04-19
3   15 1.081334 0.5405033 0.4998486 -0.9677076 2022-04-20
4   16 1.070063 0.5348910 0.4998687 -0.9729603 2022-04-21
5   17 1.072671 0.5362001 0.4998738 -0.9679795 2022-04-22
6   18 1.065505 0.5326116 0.4998677 -0.9529303 2022-04-23

画图输出R0的预测值,从预测结果看,在5月初R0是最高点,是扩散最严重的时刻,北京市也是在4月底发紧急通知,封闭朝阳区的小区,并且5月份限制办公聚集,居家上班等强制干预行动,每天进行核酸检测。到6月初才放开跨区工作,逐步恢复工作秩序。

> plot(R0~date,data=RODF)

单独看一个R0的图并不明显,我们把累计确诊,每日新增确诊和R0值,合并在一起进行输出结果, 以x轴为时间统一对齐。

> par(mfrow = c(3, 1))
> plot(cum~date,data=bj202204,type='b',col="blue")
> plot(daily~date,data=bj202204,type='b',col="red")
> 
> d2<-data.frame(date=RODF$date,R0=RODF$R0)
> d1<-data.frame(date=as.Date('2022-04-01')+0:16,R0=NA)
> d2<-rbind(d1,d2)
> plot(R0~date,data=d2,type='b',col="orange")

从2022年04月01日开始,最上面蓝色图为每日累计确诊病例,中间红色图为每日新增确诊,最下面橙色图为R0值。

通过观察R0值,很明显的能看到R0值在5月初达到顶点后,R0值开始大幅下降。从5月每日确诊人数来看趋势是平稳的,而且也在降低,因此说明北京本次防控是有效的。6月份逐步放开管制后,大家的工作和生活开始逐步恢复正常状态。

文本完整的代码,已经上传github,可以自由下载使用:https://github.com/bsspirit/infect/blob/main/code/beijing.r

转载请注明出处:
http://blog.fens.me/r-coivd19-beijing

打赏作者