我正在尝试使用Monte-Carlo模拟实现以下统计测试。此方法基于以下论文:https://journals.ametsoc.org/doi/full/10.1175/jcli4217.1
详细信息:
上文使用Monte Carlo Simulation计算了两个时期之间平均值的均值: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))
关于如何正确执行此操作的任何建议?
我将感谢有关此的任何帮助。
此代码应解决您的问题。如果您必须为大量数据进行操作,则很容易与" foreach"one_answers"多巴罗"的软件包并行。此函数获取您的数据并为数据的两个图块制作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)