我想保存此程序中的仿真结果

Question

# library (energy)
RR=100
n=10
a=2
b=4
miu1=2
miu2=4
m22=(b^2)*(1-(rho^2))
# This is the point where am having problem
# I want the programme to retain the results average0.1, average0.05 and  
# average0.01 for every 'rho' from the rho_list used for the simulation 
# but I am stuck because I don't know how to get the result
rho_list=c(0,0.3,0.6)
for (rho in rho_list){
  energy=rep(NA,RR)
  for (i in 1:RR){
    z1=rnorm(n,0,1)
    z2=rnorm(n,0,1)
    x1=miu1+a*z1
    x2=miu2+(rho*b*z1)+(sqrt(m22)*z2)
    X=matrix(c(x1,x2),byrow=TRUE,ncol=2)
    energy[i]=mvnorm.etest(X)$p.value
  }
  average0.1=sum(energy<=0.1)/RR
  average0.05=sum(energy<=0.05)/RR
  average0.01=sum(energy<=0.01)/RR
}

我希望该程序保留用于模拟的rho_list的每个rho的结果average0.1，average0.05和average0.01但是我之所以被困，是因为我不知道如何获得结果

Answer 1

您的示例不可再现，因此我给您一些模拟数据以演示如何输出结果。

rho_list=c(0,0.3,0.6)
result <- sapply(rho_list, FUN = function(rho, ...) {
    average0.1 = runif(1)
    average0.05 = runif(1)
    average0.01 = runif(1)
    c(rho = rho, a01 = average0.1, a0.05 = average0.05, a0.01 = average0.01)
}, RR = RR, n = n, a = a, b = b, miu1 = miu1, miu2 = miu2, m22 = m22, simplify = FALSE)
do.call("rbind", result)
     rho       a01      a0.05      a0.01
[1,] 0.0 0.0136175 0.08581583 0.07171591
[2,] 0.3 0.8334469 0.42103038 0.07857328
[3,] 0.6 0.8231120 0.40647485 0.65408540

Answer 2

一个选项是将结果存储在每个值的每个值的列表中，然后将它们绑定到单个数据框架中。这是一个例子。请注意，由于rho在设置代码中未定义，因此我在循环中将m22的定义替换为m22。另外，我设置了RR = 10，以节省运行代码的时间。

library(energy)
RR=10
n=10
a=2
b=4
miu1=2
miu2=4
rho_list=c(0, 0.3, 0.6)
energy_threshold = c(0.1, 0.05, 0.01) # Store energy thresholds in a vector
# Create a list of data frames. Each data frame contains the result for each
# of the three energy thresholds for one value of rho.
results = lapply(rho_list, function(rho) {
  energy=rep(NA,RR)
  for (i in 1:RR) {
    z1=rnorm(n,0,1)
    z2=rnorm(n,0,1)
    x1=miu1+a*z1
    x2=miu2+(rho*b*z1)+(sqrt((b^2)*(1-(rho^2)))*z2)
    X=matrix(c(x1,x2),byrow=TRUE,ncol=2)
    energy[i]=mvnorm.etest(X)$p.value
  }
  data.frame(rho, energy_threshold, result=sapply(energy_threshold, function(y) sum(energy <= y)/RR))
})
# Bind the three data frames into a single data frame
results = do.call(rbind, results)

这是输出：

results
  rho energy_threshold result
1 0.0             0.10    0.1
2 0.0             0.05    0.0
3 0.0             0.01    0.0
4 0.3             0.10    0.2
5 0.3             0.05    0.1
6 0.3             0.01    0.0
7 0.6             0.10    0.0
8 0.6             0.05    0.0
9 0.6             0.01    0.0

Answer 3

我将循环中的变量存储到numeric向量中，然后使用cbind()存储结果。这是整个代码：

library(energy)
RR=10
n=10
a=2
b=4
miu1=2
miu2=4
m22=(b^2)*(1-(rho^2))
average0.1 <- as.numeric()
average0.05 <- as.numeric()
average0.01 <- as.numeric()
# This is the point where am having problem
# I want the programme to retain the results average0.1, average0.05 and  
# average0.01 for every 'rho' from the rho_list used for the simulation 
# but I am stuck because I dont know how to get the result
rho_list=c(0,0.3,0.6)
for (rho in unique(rho_list)){
  energy=rep(NA,RR)
  for (i in 1:RR){
    z1=rnorm(n,0,1)
    z2=rnorm(n,0,1)
    x1=miu1+a*z1
    x2=miu2+(rho*b*z1)+(sqrt(m22)*z2)
    X=matrix(c(x1,x2),byrow=TRUE,ncol=2)
    energy[i]=mvnorm.etest(X)$p.value
}
  average0.1=rbind(average0.1, sum(energy<=0.1)/RR)
  average0.05=rbind(average0.05, sum(energy<=0.05)/RR)
  average0.01=rbind(average0.01, sum(energy<=0.01)/RR)
}