如何複製我的模擬研究

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）。