Моделирование методом Монте-Карло с использованием R: проблема с сортировкой и значимостью

Я пытаюсь реализовать следующий статистический тест, используя моделирование Монте-Карло. Этот метод основан на следующем документе: 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))

Есть предложения, как это сделать правильно?

Буду признателен за любую помощь в этом отношении.


r montecarlo significance
person Lyndz    schedule 24.05.2019    source источник
comment
Как реализовать третье условие? Какой алгоритм?   -  person Roland    schedule 24.05.2019
comment
@ Роланд .. Я тоже не уверен насчет 3-го условия. Это только что упомянуто, исходя из ранга первоначального значения разницы. Прошу прощения за непонимание относительно этого теста. Я также впервые делаю тест Монте-Карло.   -  person Lyndz    schedule 24.05.2019
comment
Я уже пытался помочь вам: задайте конкретный вопрос на перекрестной проверке, чтобы понять это предложение. Тогда возвращайтесь сюда, если вам понадобится помощь с Р.   -  person Roland    schedule 24.05.2019

Ответы (1)


arrow_upward
1
arrow_downward

Этот код должен решить вашу проблему. Если вам нужно сделать это для большого количества данных, это легко распараллелить с помощью пакетов «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)
person Santiago I. Hurtado    schedule 24.05.2019
comment
Большое спасибо! - person Lyndz; 24.05.2019