Как повторить мое исследование симуляции

Question

Я новичок в R, и я занимаюсь симуляцией, и мне удалось создать один образец, который я хочу сделать. Однако я понятия не имею, как мне повторить то, что я сделал.

Вот программа, которую я написал до сих пор:

I <- 500       # number of observations
J <- 18        # total number of items
K <- 6         # number of testlets
JK <-3         # number of items within a testlet
response <- matrix(0, I, J)  # null binary (0, 1) response matrix 
unit <- matrix(1, JK, 1)     # unit vector

set.seed(1234)

# Multidimensional 3-pl model
pij <- function(a,b,c,theta,gamma) {c+(1-c)*(1/(1+exp(-1.7*a*(theta-b-gamma))))}

# Assigning a and b parameter values
a <- c(.8,.9,.7,.8,.9,.7,.8,.9,.7,.8,.9,.7,.8,.9,.7,.8,.9,.7)
b <-c(1,0,-1.5,1,0,-1.5,1,0,-1.5,1,0,-1.5,1,0,-1.5,1,0,-1.5)
# Assigning c-parameter, each 3 items (c-parameter & testlet effect)
#(small&small, small&large, large&small, large&large, mixed&small, mixed&large)
c <- c(.2,.2,.2,.2,.2,.2,.5,.5,.5,.5,.5,.5,.2,.33,.5,.2,.33,.5)    

theta <- rnorm(I, 0, 1)   # random sampling theta-values from normal dist. M=0, SD=1

gamma1 <- rnorm(I, 0, .2)  # small testlet effect: random sampling gamma from normal dist. M=0, SD=.2
gamma2 <- rnorm(I, 0, 1)   # large testlet effect: random sampling gamma from normal dist. M=0, SD=1
gamma3 <- rnorm(I, 0, .2)  # small testlet effect: random sampling gamma from normal dist. M=0, SD=.2
gamma4 <- rnorm(I, 0, 1)   # large testlet effect: random sampling gamma from normal dist. M=0, SD=1
gamma5 <- rnorm(I, 0, .2)  # small testlet effect: random sampling gamma from normal dist. M=0, SD=.2
gamma6 <- rnorm(I, 0, 1)   # large testlet effect: random sampling gamma from normal dist. M=0, SD=1

# implementing that the testlet effect is same for the items within a testlet
gamma1T <- gamma1 %*% t(unit)
gamma2T <- gamma2 %*% t(unit)
gamma3T <- gamma3 %*% t(unit)
gamma4T <- gamma4 %*% t(unit)
gamma5T <- gamma5 %*% t(unit)
gamma6T <- gamma6 %*% t(unit)

gammaT <- matrix(c(gamma1T, gamma2T, gamma3T, gamma4T, gamma5T, gamma6T), I, J)  # getting all the gammas together in a large matrix

# Generating data using the information above
for(i in 1:I) {
  for(j in 1:J) {
    response[i, j] <- ifelse(pij(a=a[j], b=b[j], c=c[j], theta=theta[i], gamma=gammaT[i,j]) < runif(1), 0, 1)
  }
}

Таким образом, я получаю набор данных «ответ». Я хочу повторить это и получить, скажем, 1000 «ответных» наборов данных. Я думаю, что это можно сделать путем репликации случайной выборки для «тета» и «гамма», но у меня нет идеи сделать это на самом деле.

Большое, большое спасибо заранее,

Hanjoe.

Answer 1

Совет Stedy разумен, за исключением одного: НЕ увеличивать начальное число в цикле for.

Как я понимаю предложение Стеди, set.seed(i) будет вызываться в цикле for для каждой симуляции, при этом i увеличивается на каждой итерации. Эта стратегия известна для введения (потенциально большого) смещения из-за корреляций между сгенерированными последовательностями.

Вместо этого устанавливает начальное значение один раз в начале, то есть перед циклом for. Например, вы можете использовать текущее время Unix в качестве начального числа или прочитать его из файла со случайными числами (например, из random.org ). Кроме того, не забудьте сохранить семена с вашими результатами, например, распечатать его в файл журнала. Если вы хотите снова воспроизвести точные результаты предыдущего набора репликаций, вам просто нужно установить соответствующее начальное число.

Если вы хотите, чтобы другие могли точно копировать ваши результаты, вам также следует указать версию R (и, возможно, операционную систему), которую вы использовали (поскольку реализации RNG могут различаться).

Кроме того, имитационная репликация является параллельной задачей , то есть вы можете легко выполнять репликацию параллельно, если у вас есть многоядерный компьютер (например, с rparallel ). В этом случае, однако, необходима особая осторожность в отношении случайных чисел (например, см. этот документ ).

Answer 2

Я бы взял локальные переменные и превратил их в функцию. Затем создайте цикл for(), вызывайте функцию и увеличивайте set.seed() на единицу каждый раз, когда вызывается функция для длины цикла for().