如何复制我的模拟研究

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个“响应”数据集。 我认为这可以通过复制“theta”和“gamma”的随机采样来完成，但是我没有想法实际做到这一点。

很多很多预先感谢，

Hanjoe。

Answer 1

我会采取局部变量，并把它们变成一个函数。然后创建一个for()循环，调用该函数，并在每次调用函数for()循环的长度时递增set.seed()。

Answer 2

Stedy的建议是合理的，除了一件事：不要增加for循环中的种子。

就我所了解的Stedy的建议而言，set.seed(i)将在for循环内对每个模拟进行调用，每个迭代中增加i。该策略是known，以引入由于所生成的序列之间的相关性而可能较大的偏差。

而是，在开始时即在for循环之前设置种子一次。例如，您可以使用当前的Unix时间作为种子，或者使用随机数从文件中读取一个（例如，从random.org）。另外，请确保将结果存储在种子中，例如将其打印到日志文件。如果您想再次复制上一组复制的准确结果，您只需设置相应的种子。

如果您希望其他人能够完全复制您的结果，您还应该指定您使用的R版本（也可能是操作系统）（因为RNG实施可能会有所不同）。另外，仿真复制是embarassingly parallel任务，即如果您有多核机器（例如rparallel），则可以并行轻松地执行复制。但是，在这种情况下，需要特别注意随机数（例如，参见this paper）。