使用 R 的蒙特卡洛模拟：排序和重要性问题答案

【问题标题】：Monte Carlo Simulation using R: Problem with sorting and significance使用 R 的蒙特卡洛模拟：排序和重要性问题
【发布时间】：2019-10-10 18:46:58
【问题描述】：

我正在尝试使用蒙特卡罗模拟来实现以下统计测试。该方法基于以下论文： https://journals.ametsoc.org/doi/full/10.1175/JCLI4217.1

详情：

上述论文使用蒙特卡罗模拟计算了两个时期：1961-1983 和 1984-2000 热带气旋通过频率（非正态分布）的均值差异的显着性。

这应该是一个双尾测试。

提供了以下步骤：

1)。首先，准备了9999个随机排序的40年台风通过频率时间序列。

2)。计算前 23 年值（1961-1983）的平均值减去后 17 年值的平均值。

3)。从10000个样本中原始差值的等级来估计显着性水平。

这是我目前所拥有的

假设我有以下数据集。列表示每年的计数，而行表示经纬度坐标（为简单起见，此处为数字）。

A<-matrix(floor(runif(100,min=0,max=20)),nrow=5,ncol=40)
colnames(A)<-c("X1961","X1962","X1963","X1964","X1965","X1966","X1967","X1968","X1969","X1970","X1971","X1972","X1973","X1974","X1975","X1976","X1977","X1978","X1979","X1980","X1981","X1982","X1983","X1984","X1985","X1986","X1987","X1988","X1989","X1990","X1991","X1992","X1993","X1994","X1995","X1996","X1997","X1998","X1999","X2000")

set.seed(1)
rand <- sample(nrow(A),9999,replace=TRUE)
A[rand,]

问题（更新）

我对如何在 R 中正确执行此操作感到困惑。我应该每行执行蒙特卡罗测试。所以在一行中这样做：

A[rand[1],]

X1961 X1962 X1963 X1964 X1965 X1966 X1967 X1968 X1969 X1970 X1971 X1972 
X1973 
5    14    11    17    16    17    11     2     8     3    13    10     
1 
X1974 X1975 X1976 X1977 X1978 X1979 X1980 X1981 X1982 X1983 X1984 X1985 
X1986 
10    15     5     3     6    15    19     5    14    11    17    16    
17 
X1987 X1988 X1989 X1990 X1991 X1992 X1993 X1994 X1995 X1996 X1997 X1998 
X1999 
11     2     8     3    13    10     1    10    15     5     3     6    
15 
X2000 
19

原文：

A[1,]

X1961 X1962 X1963 X1964 X1965 X1966 X1967 X1968 X1969 X1970 X1971 X1972 
X1973 
18     1     6     7     3    12    19     0    17    17     0    10    
16 
X1974 X1975 X1976 X1977 X1978 X1979 X1980 X1981 X1982 X1983 X1984 X1985 
X1986 
3     4     0    15     8    17     1    18     1     6     7     3    
12 
X1987 X1988 X1989 X1990 X1991 X1992 X1993 X1994 X1995 X1996 X1997 X1998 
X1999 
19     0    17    17     0    10    16     3     4     0    15     8    
17 
X2000 
1

预期输出*

我想在这个测试的原始矩阵中添加一个 pvalue 列。显着性检验应按行进行。当然，这可以通过使用 apply() 函数来实现。

问题

如何实现第三个条件？另外，蒙特卡洛测试中步骤 1 的顺序是否重要？

我觉得我误解了第 1 步，我应该为此使用 replicate() 吗？像这样？

rand<-replicate(40,sample(nrow(A),9999,replace=T))

关于如何正确执行此操作的任何建议？

对于这方面的任何帮助，我将不胜感激。

【问题讨论】：

"如何实现第三个条件？"算法是什么？
@Roland..我也不确定第三种情况。它刚刚提到，“从原始差值的等级”。对于对这次测试缺乏了解，我深表歉意。这也是我第一次做蒙特卡罗测试。
我已经尝试指导您：在 Cross Validated 上提出一个特定问题以理解这句话。如果您需要 R 方面的帮助，请回到这里。

标签： r montecarlo significance

【解决方案1】：

此代码应该可以解决您的问题。如果您必须为大量数据执行此操作，它很容易与包“foreach”和“doParallel”并行化。此函数获取您的数据并为数据的两个图块制作 nrep 样本，然后获取均值的差异。以此计算均值差的FDP，然后查看均值数据差的百分位数得到p值。

  my.fun <- function(x,nrep = 1000,breakpoint){
    # x is the data
    # nrep is the amount of simulations
    # breakpoint is where the breakpoint is
    set.seed(12345)
    a_sim <- vector(mode = 'double', length = nrep)
    n <- length(x)
    for(i in 1:nrep){
      aux1 <- sample(x,size = breakpoint,replace = T)
      aux2 <- sample(x,size = n-breakpoint,replace = T)
      a_sim[i] <- abs(mean(aux1) - mean(aux2))
    }
    cum_dist_func <- ecdf(a_sim)
    p <- 1-cum_dist_func(abs(mean(x[1:breakpoint])-mean(x[(breakpoint+1):n])))

    return(p)
  }


  pvalue <- apply(X = A,MARGIN = 1,FUN = my.fun,breakpoint = 23 )
  A <- cbind(A,pvalue)

【讨论】：

非常感谢！