R: 使用Mersenne-Twister种子的runif生成的随机值严重聚集

Question

R: 使用Mersenne-Twister种子的runif生成的随机值严重聚集

5

我们在使用R的runif函数并使用set.seed命令时，遇到了一个奇怪的情况，这个函数使用了kind=NULL参数设置了种子（如果我没有理解错误的话，它等效于kind="default"，默认的是"Mersenne-Twister"）。

在调用runif函数之前，我们使用由上游系统生成的(8位数)唯一ID来设置种子。

seeds = c(
  "86548915", "86551615", "86566163", "86577411", "86584144", 
  "86584272", "86620568", "86724613", "86756002", "86768593", "86772411", 
  "86781516", "86794389", "86805854", "86814600", "86835092", "86874179", 
  "86876466", "86901193", "86987847", "86988080")

random_values = sapply(seeds, function(x) {
  set.seed(x)
  y = runif(1, 17, 26)
  return(y)
})

这会给出非常集中的值。

> summary(random_values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  25.13   25.36   25.66   25.58   25.83   25.94

当我们使用kind = "Knuth-TAOCP-2002"时，runif的这种行为就会消失，我们得到的值看起来更加均匀分布。

random_values = sapply(seeds, function(x) {
  set.seed(x, kind = "Knuth-TAOCP-2002")
  y = runif(1, 17, 26)
  return(y)
})

输出被省略。

最有趣的是，这种情况只会在Ubuntu上发生，而不会在Windows上发生（下面是Ubuntu和Windows的sessionInfo输出）。

Windows输出：

> seeds = c(
+   "86548915", "86551615", "86566163", "86577411", "86584144", 
+   "86584272", "86620568", "86724613", "86756002", "86768593", "86772411", 
+   "86781516", "86794389", "86805854", "86814600", "86835092", "86874179", 
+   "86876466", "86901193", "86987847", "86988080")
> 
> random_values = sapply(seeds, function(x) {
+   set.seed(x)
+   y = runif(1, 17, 26)
+   return(y)
+ })
> 
> summary(random_values)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  17.32   20.14   23.00   22.17   24.07   25.90

有人能帮忙理解发生了什么吗？

Ubuntu

R version 3.4.0 (2017-04-21)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 16.04.2 LTS

Matrix products: default
BLAS: /usr/lib/libblas/libblas.so.3.6.0
LAPACK: /usr/lib/lapack/liblapack.so.3.6.0

locale:
[1] LC_CTYPE=en_US.UTF-8          LC_NUMERIC=C                 
 [3] LC_TIME=en_US.UTF-8           LC_COLLATE=en_US.UTF-8       
 [5] LC_MONETARY=en_US.UTF-8       LC_MESSAGES=en_US.UTF-8      
 [7] LC_PAPER=en_US.UTF-8          LC_NAME=en_US.UTF-8          
 [9] LC_ADDRESS=en_US.UTF-8        LC_TELEPHONE=en_US.UTF-8     
[11] LC_MEASUREMENT=en_US.UTF-8    LC_IDENTIFICATION=en_US.UTF-8

attached base packages:
[1] parallel  stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
[1] RMySQL_0.10.8               DBI_0.6-1                  
 [3] jsonlite_1.4                tidyjson_0.2.2             
 [5] optiRum_0.37.3              lubridate_1.6.0            
 [7] httr_1.2.1                  gdata_2.18.0               
 [9] XLConnect_0.2-12            XLConnectJars_0.2-12       
[11] data.table_1.10.4           stringr_1.2.0              
[13] readxl_1.0.0                xlsx_0.5.7                 
[15] xlsxjars_0.6.1              rJava_0.9-8                
[17] sqldf_0.4-10                RSQLite_1.1-2              
[19] gsubfn_0.6-6                proto_1.0.0                
[21] dplyr_0.5.0                 purrr_0.2.4                
[23] readr_1.1.1                 tidyr_0.6.3                
[25] tibble_1.3.0                tidyverse_1.1.1            
[27] rBayesianOptimization_1.1.0 xgboost_0.6-4              
[29] MLmetrics_1.1.1             caret_6.0-76               
[31] ROCR_1.0-7                  gplots_3.0.1               
[33] effects_3.1-2               pROC_1.10.0                
[35] pscl_1.4.9                  lattice_0.20-35            
[37] MASS_7.3-47                 ggplot2_2.2.1              

loaded via a namespace (and not attached):
[1] splines_3.4.0      foreach_1.4.3      AUC_0.3.0          modelr_0.1.0      
 [5] gtools_3.5.0       assertthat_0.2.0   stats4_3.4.0       cellranger_1.1.0  
 [9] quantreg_5.33      chron_2.3-50       digest_0.6.10      rvest_0.3.2       
[13] minqa_1.2.4        colorspace_1.3-2   Matrix_1.2-10      plyr_1.8.4        
[17] psych_1.7.3.21     XML_3.98-1.7       broom_0.4.2        SparseM_1.77      
[21] haven_1.0.0        scales_0.4.1       lme4_1.1-13        MatrixModels_0.4-1
[25] mgcv_1.8-17        car_2.1-5          nnet_7.3-12        lazyeval_0.2.0    
[29] pbkrtest_0.4-7     mnormt_1.5-5       magrittr_1.5       memoise_1.0.0     
[33] nlme_3.1-131       forcats_0.2.0      xml2_1.1.1         foreign_0.8-69    
[37] tools_3.4.0        hms_0.3            munsell_0.4.3      compiler_3.4.0    
[41] caTools_1.17.1     rlang_0.1.1        grid_3.4.0         nloptr_1.0.4      
[45] iterators_1.0.8    bitops_1.0-6       tcltk_3.4.0        gtable_0.2.0      
[49] ModelMetrics_1.1.0 codetools_0.2-15   reshape2_1.4.2     R6_2.2.0          
[53] knitr_1.15.1       KernSmooth_2.23-15 stringi_1.1.5      Rcpp_0.12.11

Windows

> sessionInfo()
R version 3.3.2 (2016-10-31)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows >= 8 x64 (build 9200)

locale:
[1] LC_COLLATE=English_India.1252  LC_CTYPE=English_India.1252    LC_MONETARY=English_India.1252
[4] LC_NUMERIC=C                   LC_TIME=English_India.1252    

attached base packages:
[1] graphics  grDevices utils     datasets  grid      stats     methods   base     

other attached packages:
 [1] bindrcpp_0.2         h2o_3.14.0.3         ggrepel_0.6.5        eulerr_1.1.0         VennDiagram_1.6.17  
 [6] futile.logger_1.4.3  scales_0.4.1         FinCal_0.6.3         xml2_1.0.0           httr_1.3.0          
[11] wesanderson_0.3.2    wordcloud_2.5        RColorBrewer_1.1-2   htmltools_0.3.6      urltools_1.6.0      
[16] timevis_0.4          dtplyr_0.0.1         magrittr_1.5         shiny_1.0.5          RODBC_1.3-14        
[21] zoo_1.8-0            sqldf_0.4-10         RSQLite_1.1-2        gsubfn_0.6-6         proto_1.0.0         
[26] gdata_2.17.0         stringr_1.2.0        XLConnect_0.2-12     XLConnectJars_0.2-12 data.table_1.10.4   
[31] xlsx_0.5.7           xlsxjars_0.6.1       rJava_0.9-8          readxl_0.1.1         googlesheets_0.2.1  
[36] jsonlite_1.5         tidyjson_0.2.1       RMySQL_0.10.9        RPostgreSQL_0.4-1    DBI_0.5-1           
[41] dplyr_0.7.2          purrr_0.2.3          readr_1.1.1          tidyr_0.7.0          tibble_1.3.3        
[46] ggplot2_2.2.0        tidyverse_1.0.0      lubridate_1.6.0     

loaded via a namespace (and not attached):
 [1] gtools_3.5.0         assertthat_0.2.0     triebeard_0.3.0      cellranger_1.1.0     yaml_2.1.14         
 [6] slam_0.1-40          lattice_0.20-34      glue_1.1.1           chron_2.3-48         digest_0.6.12.1     
[11] colorspace_1.3-1     httpuv_1.3.5         plyr_1.8.4           pkgconfig_2.0.1      xtable_1.8-2        
[16] lazyeval_0.2.0       mime_0.5             memoise_1.0.0        tools_3.3.2          hms_0.3             
[21] munsell_0.4.3        lambda.r_1.1.9       rlang_0.1.1          RCurl_1.95-4.8       labeling_0.3        
[26] bitops_1.0-6         tcltk_3.3.2          gtable_0.2.0         reshape2_1.4.2       R6_2.2.0            
[31] bindr_0.1            futile.options_1.0.0 stringi_1.1.2        Rcpp_0.12.12.1

- tchakravarty

2

“上游过程”是生成种子的过程。这似乎是由您提供的特定一组种子的产物，并不适用于类似规模的随机种子集合。 - dww

如果我在Windows上运行您的代码，我会得到相同的结果。我在Windows 10上使用了3.4.1。通常情况下，不同平台之间的随机数生成没有区别。除非我误解了您所说的“这种情况不会发生”的确切含义。 - MrFlick

@MrFlick，我可以提供你在我的Windows 10机器上运行R 3.3.2的sessionInfo，在这台机器上随机数并没有像这个问题中那样呈非随机分布的情况。 - tchakravarty

此外，这被踩的有点奇怪。 - tchakravarty

@tchakravarty 是的。提供会话信息，以便在Windows上显示不同的结果，并展示这些结果。使用相同的种子，所有版本的R应该能够再现结果。 - MrFlick

显示剩余5条评论

3个回答

2

作为证明您的序列是边缘情况的更多证据，您可以关注构建的看似随机值的范围。17和26有点分散注意力。在0和1上使用均匀分布重复您的实验会产生同样不太可能的结果。

f <- function(x) {
  set.seed(x)
  runif(1)
}

 check_range <-function(seeds){
   vals <- sapply(seeds,f)
   max(vals)-min(vals)
}

当对您的种子运行时：

> check_range(seeds)
[1] 0.09026112

在运行21个随机种子时，对check_range(seeds)进行合理的建模是它是从U(0,1)中抽取大小为21的随机样本的样本范围。它的理论密度如下：

f <- function(x){420*x^19*(1-x)}

我们可以使用这个来计算观察到0.09或更小范围的概率：

> integrate(f,0,0.09)
2.334272e-20 with absolute error < 2.6e-34

为了验证在使用Mersenne Twister作为随机数生成器时，用样本范围作为种子是否合理，你可以进行以下实验：

ranges <- replicate(1000,check_range(sample(8548915:86988080,21)))
x <- seq(0,1,0.01)
y <- f(x)
hist(ranges,freq = FALSE,xlim =c(0,1))
points(x,y,type = "l")
abline(v=0.09)

输出:

密度直方图与理论密度相当吻合。你问题中的21个种子集合代表了一个极端的离群值。这很不可能是巧合造成的，也不太可能是Mersenne Twister中某些潜在缺陷所导致的。最有可能的解释是Mersenne Twister本身参与了产生这些21个值的过程（但当然不是简单地使用sample()来绘制21个值）。

- John Coleman

1

当您使用单个种子的Mersenne Twister时，可以合理地假设生成的值大致上是独立且同分布的。不幸的是，无法保证从两个以不同种子开始的流中生成的值。例如，请参见this SC thread。

在您的情况下，我建议使用SC线程中建议的其中一种种子选择策略，或者切换到具有更好并行流保证的PRNG。一个选项是L'Ecuyer的“RngStreams”生成器：

set.seed(0, kind = "L'Ecuyer-CMRG")

即使使用了那个伪随机数生成器，我也不确定是否仍然可以使用任意种子来初始化PRNG并获得大致独立的流。

至于Ubuntu和Windows之间的差异，有可能其中一个系统使用32位生成器，而另一个使用64位。

- Patrick Perry

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Ista · Accepted Answer

注意：本回答总结了在R-devel邮件列表上讨论此问题的要点。我只是试图捕捉和总结最初在那里表达的想法。

尽管您保证这些数字不是特别构造的边缘情况，但它们似乎确实是这样。以下是原始序列以及检查生成值分布的代码：

seeds = c(
    86548915, 86551615, 86566163, 86577411, 86584144, 86584272,
    86620568, 86724613, 86756002, 86768593, 86772411, 86781516,
    86794389, 86805854, 86814600, 86835092, 86874179, 86876466,
    86901193, 86987847, 86988080)
checkit <- function(seeds) {
    sapply(seeds, function(x) {
        set.seed(x)
        y = runif(1, 17, 26)
        return(y)
    })}

原始序列变异很小，正如所指出的：

  summary(checkit(seeds+0))
  ## Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  ##25.13   25.36   25.66   25.58   25.83   25.94

原始序列似乎具有某种特殊性质，因为对它进行最小限度的修改并不会产生相同令人惊讶的结果：

summary(checkit(seeds+1))
## Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
## 17.18   19.65   22.75   22.02   24.37   25.79

summary(checkit(seeds-1))
## Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##17.15   18.44   19.92   20.77   22.97   25.95

在原始序列所覆盖的所有种子中，预期数量会产生观察范围内的值：

possible.seeds <- min(seeds):max(seeds)

s25 <- Filter(function(s){
    set.seed(s)
    x <- runif(1,17,26)
    x > 25.12 & x < 25.95},
    possible.seeds)

length(s25)/length(possible.seeds)
##[1] 0.09175801

然而，原始序列中的所有值都在这个子集中（当然，我们已经知道这一点...）。

table(seeds %in% s25)

##TRUE 
##  21

所有这些都指向一个可能性，即原始序列实际上是一个（也许是无意中）特别构造的边缘情况。