131-R茶话会23-R的随机数有点坑

北野茶缸子

发布于 2022-05-19 11:34:13

66300

代码可运行

文章被收录于专栏：北野茶缸子的专栏北野茶缸子的专栏

运行总次数：0

代码可运行

前言

最近我在复现一篇文章的操作。发现每一次生成的结果都有所不同。

难道是我的操作出了问题？难道是我用的R 包版本不对，函数不同？难道是随机数的问题？

后来发现，果然是随机数的问题。记得之前[[103-R茶话会18-随机数和取子集是天生不和吗？]] 就曾经聊过。

1-随机数，老是变

通常我们都会通过set.seed 来设置随机数。

R 内内置了许多的随机相关的函数，比如：

> runif(1)
[1] 0.06532152
> runif(2)
[1] 0.2081815 0.8665349
> sample(10,3)
[1]  4  9 10
> rnorm(3)
[1]  0.18983511  0.09097173 -0.10174577
> rnorm(3)
[1]  0.89563748 -0.01084824  2.07215441

这里就不详细展开介绍了。

电脑永远无法真正实现随机数，R 也不例外，关于生成看似“随机”数的原理，可以参考：Set.seed in R - Control Random Numbers - ProgrammingR[1]

如果我们想要控制随机数，或者说，其他人可以重复我们执行涉及到随机数的函数，可以用内置的set.seed 提前声明：

> set.seed(33)
> runif(2)
[1] 0.4459405 0.3946503
> set.seed(33)
> runif(2)
[1] 0.4459405 0.3946503

所以说，一个set.seed，我们就可以一劳永逸的不管不顾了吗？

2-随机数，即用即换

我们如果希望随机函数生成指定结果，永远要在其之前配置相关的种子。

> set.seed(33)
> runif(3)
[1] 0.4459405 0.3946503 0.4837289
> runif(3)
[1] 0.9188760 0.8438814 0.5173496
> set.seed(33)
> runif(3)
[1] 0.4459405 0.3946503 0.4837289

换句话说，一旦某个随机函数触发，设置的种子就会重置。

即使更换其他函数，也不例外：

> set.seed(33)
> runif(3)
[1] 0.4459405 0.3946503 0.4837289
> set.seed(33)
> runif(3)
[1] 0.4459405 0.3946503 0.4837289
> set.seed(33)
> rnorm(4)
[1] -0.13592452 -0.04079697  1.01053901 -0.15826244
> runif(3)
[1] 0.01551696 0.11799116 0.69098590

这实在有点麻烦啊。

如果是执行两次或多次set.seed，可以累积计算吗？

> set.seed(33)
> set.seed(33)
> runif(3)
[1] 0.4459405 0.3946503 0.4837289
> runif(3)
[1] 0.9188760 0.8438814 0.5173496

显然不行。

换句话说，如果你的循环调用了随机函数，则这个seed 一定要设定在循环内部，否则你永远无法重复出自己的代码。

就比如我遇到的项目：

sce = lapply(unique(batchFactor) , function(current.batch){
  idx = which(batchFactor == current.batch)
  current.sce = SingleCellExperiment(assays = list("counts" = counts[,idx]), colData = meta[idx,])
  clusters <- quickCluster(current.sce, method="igraph", use.ranks=TRUE, d=50, min.mean=0.1, min.size = 50)
  current.sce <- computeSumFactors(current.sce, clusters=clusters)
  return(current.sce)
})
sce = do.call(multiBatchNorm , sce)
sce = do.call(cbind, sce)

最终发现，原来这段代码调用了随机函数。难怪我重复不出来！

当我在循环内部加了一个种子：

sce = lapply(unique(batchFactor) , function(current.batch){
  set.seed(32)
  idx = which(batchFactor == current.batch)
  current.sce = SingleCellExperiment(assays = list("counts" = counts[,idx]), colData = meta[idx,])
  clusters <- quickCluster(current.sce, method="igraph", use.ranks=TRUE, d=50, min.mean=0.1, min.size = 50)
  current.sce <- computeSumFactors(current.sce, clusters=clusters)
  return(current.sce)
})
sce = do.call(multiBatchNorm , sce)
sce = do.call(cbind, sce)

> sce2class: SingleCellExperiment 
dim: 10000 4645 
metadata(0):
assays(2): counts logcounts
rownames(10000): RPL21 TMSB4X ... MIR637 PRAMEF13
rowData names(0):
colnames(4645): Cy72_CD45_H02_S758_comb Cy72_CD45_D09_S717_comb ...
  CY94_CD45NEG_CD90POS_2_F04_S64_comb
  CY94_CD45NEG_CD90POS_2_D08_S44_comb
colData names(5): cell sample malignancy celltype sizeFactor
reducedDimNames(0):
altExpNames(0):

虽然一般来说，我们并不会在循环内部使用相同的种子。但比如许多的图聚类算法，都会涉及到随机。那么如果不在lapply 中特别声明，或者仅仅在全局声明了一次set.seed，那么你就等着无法重复的随机结果吧：

> set.seed(32)
> lapply(1:3, function(x) runif(3))
[[1]]
[1] 0.5058405 0.5948084 0.8087471

[[2]]
[1] 0.7288197 0.1519876 0.9561873

[[3]]
[1] 0.7535377 0.8520623 0.6734418

你必须得在循环内部声明：

> lapply(1:3, function(x) {set.seed(32); runif(3) })
[[1]]
[1] 0.5058405 0.5948084 0.8087471

[[2]]
[1] 0.5058405 0.5948084 0.8087471

[[3]]
[1] 0.5058405 0.5948084 0.8087471

3-随机数的触发与设定都是全局的

关于随机数的触发，以Y叔的这个文章为例：ggplot2的一个坑[2]

上面我们说了循环内部，但这种随机数的触发，还是在全局的作用域。如果是在函数中呢？

> set.seed(32)
> test <- function(x) {invisible(sample(2))}
> test(2)
> runif(3)
[1] 0.8087471 0.7288197 0.1519876
> set.seed(32)
> runif(3)
[1] 0.5058405 0.5948084 0.8087471

这也是为什么，老版本的ggplot2 加载包，会把你的随机数吃了。

同样的随机数的设置，也会默认作用在全局：

> set.seed(32)
> runif(3)
[1] 0.5058405 0.5948084 0.8087471
> runif(3)
[1] 0.7288197 0.1519876 0.9561873
> test <- function(x) {set.seed(32)}
> test(2)
> runif(3)
[1] 0.5058405 0.5948084 0.8087471

4-我似乎没有找到什么好的方法

参见：Random seed — with_seed • withr[3]

讲真话，遇到一个随机函数，就设置一个随机数，实在是痛苦。

虽然这个函数，可以显示的在函数内部调用随机数和函数，帮助我们明确随机数与随机函数的对应关系：

> withr::with_seed(32, runif(1:5)) 
[1] 0.5058405 0.5948084 0.8087471 0.7288197 0.1519876

但你依然没法帮助我解决循环啊：

> withr::with_seed(32, lapply(1:3, function(x){runif(1:5)}))
[[1]]
[1] 0.5058405 0.5948084 0.8087471 0.7288197 0.1519876

[[2]]
[1] 0.9561873 0.7535377 0.8520623 0.6734418 0.3871255

[[3]]
[1] 0.6580025 0.3213696 0.6120750 0.7726361 0.2662122

这个包里的，with_tempfile 函数我还挺喜欢的：

with_tempfile("file1", {
  print(file1)
  writeLines("foo", file1)
  readLines(file1)
})
#> [1] "/var/folders/dt/r5s12t392tb5sk181j3gs4zw0000gn/T//RtmpC6vpTp/file15b492f945411"
#> [1] "foo"

即用即丢，不污染环境。

此外，我也推荐大家在写涉及到随机数的代码时，使用withr::with_seed，起码告诉并且提醒自己，哪个种子，对应哪个随机函数。防止被吞，或者吞了其他人的种子。

额外补充

并不是所有R 的使用者都会注意到随机数的问题了。因此，考虑是否发生随机数的一个方法是，既然set.seed是即用即废，那么这也是个用来判断某个步骤是否调用了随机数的方法。

如果测试的随机数结果发生了变化，自然就可以定位到某句确实调用了随机函数：

> set.seed(32)
> runif(2)
[1] 0.5058405 0.5948084
> set.seed(32)
> runif(2)
[1] 0.5058405 0.5948084
> set.seed(32)
> sample(10)
 [1]  6  4  1  8  7  3  5  2 10  9
> runif(2)
[1] 0.2662122 0.7488021
> set.seed(32)
> print("Hello")
[1] "Hello"
> runif(2)
[1] 0.5058405 0.5948084

但这种痛苦的代码检查过程，比找bug 还难受啊！

你们有没有什么解决随机数的方案呢？

参考资料

[1]

Set.seed in R - Control Random Numbers - ProgrammingR: https://www.programmingr.com/examples/neat-tricks/sample-r-function/set-seed-in-r-control-random-numbers/

[2]

ggplot2的一个坑: https://guangchuangyu.github.io/cn/2018/01/ggplot2-seed/

[3]

Random seed — with_seed • withr: https://withr.r-lib.org/reference/with_seed.html

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-04-26，如有侵权请联系 cloudcommunity@tencent.com 删除

https

网络安全

本文分享自北野茶缸子微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

https

网络安全

登录后参与评论

0 条评论

热度