How to do simulation in my code and get the average value?

[Vilter]

library(pcalg)
## define parameters
p <- 10 # number of random variables
n <- 1000 # number of samples
szero <- 10  # expected number of edges in DAG
s <- szero / (p*(p-1)) # sparsness of the graph



## generate random data
# set.seed(42)
g <- randomDAG(p,s) # generate a random DAG
d <- rmvDAG(n,g, errDist="normal") # generate random samples

suppressWarnings({
gCPDAG <-
  pcAlgo(d,alpha=0.05, directed=TRUE) # estimate of the CPDAG
})

(shd.val <- shd(g,gCPDAG))

I want to do such above codes 12 times. That is 12 runs simulation. Then I want to get the average value of shd.val.

I also want to compute the running time of these 12 simulations.

How to do these?

[technocrat]

library(pcalg)
## define parameters
p <- 10 # number of random variables
n <- 1000 # number of samples
szero <- 10  # expected number of edges in DAG
s <- szero / (p*(p-1)) # sparsness of the graph
## generate random data
# set.seed(42)
g <- randomDAG(p,s) # generate a random DAG
set.seed(NULL) 

take_samples <- function() {
  d = rmvDAG(n,g, errDist="normal") # generate random samples
  suppressWarnings({
    gCPDAG =
      pcAlgo(d,alpha=0.05, directed=TRUE) # estimate of the CPDAG
    })
}

library(tictoc) # simple timer 
# receiver object
results <- vector(length = 12)
tic()
for(i in 1:12) results[i] = shd(g,take_samples())
results
#>  [1] 2 2 2 4 2 2 2 5 3 3 3 4
mean(results)
#> [1] 2.833333
toc()
#> 0.172 sec elapsed

^{Created on 2023-09-04 with reprex v2.0.2}

[Vilter]

How to set the seed in your solution to make the result reproducible? Thank you so much!

[technocrat]

@technocrat

rmvDAG

You could set.seed() within the take_samples() function, which would use a different seed from the value of 42 used for randomDAG() or you could remove set.seed(NULL), but either would make results contain identical values.