R Under development (unstable) (2024-02-05 r85863 ucrt) -- "Unsuffered Consequences" Copyright (C) 2024 The R Foundation for Statistical Computing Platform: x86_64-w64-mingw32/x64 R is free software and comes with ABSOLUTELY NO WARRANTY. You are welcome to redistribute it under certain conditions. Type 'license()' or 'licence()' for distribution details. R is a collaborative project with many contributors. Type 'contributors()' for more information and 'citation()' on how to cite R or R packages in publications. Type 'demo()' for some demos, 'help()' for on-line help, or 'help.start()' for an HTML browser interface to help. Type 'q()' to quit R. > library(pcalg) > suppressWarnings(RNGversion("3.5.0")) > > p <- 10 > ip <- seq_len(p) # 1:p = 1 2 .. p > > all.equal.1 <- function(x1, xx, tol = 1e-13) { # F38 Lnx works w/ 1e-15; M1 Mac + Accelerate failed w/ 1e-14 + stopifnot(length(x1) == 1L, length(xx) >= 1L) + ## the one xx[] that is close to x1 : + x. <- xx[which.min(abs(xx - x1))] + ## r528{mmaechler}: (if(abs(x1) < 1e-100) abs(x.- x1) else abs(x./x1 - 1)) <= tol + all.equal(x1, x., tolerance = tol, scale = if(abs(x1) > 0) abs(x1)) + } > > for (i in 1:50) { + cat("i=",i,"\n--==\n",sep="") # need to see where it fails (if ..) + set.seed(i) + ## generate and draw random DAG : + myDAG <- randomDAG(p, prob = 0.4) + myCPDAG <- dag2cpdag(myDAG) + mcov <- trueCov(myDAG) + + x <- sample(ip, 1) + y1 <- sample(setdiff(ip, x), 1) + y2 <- sample(setdiff(ip,c(x,y1)), 1) + y3 <- sample(setdiff(ip,c(x,y1,y2)),1) + ## plot(myCPDAG) + eff.true1 <- causalEffect(myDAG, y1, x) + eff.true2 <- causalEffect(myDAG, y2, x) + eff.true3 <- causalEffect(myDAG, y3, x) + ## cat("x=",x," y1=",y1," eff=",eff.true1,"\n") + ## cat("x=",x," y1=",y2," eff=",eff.true2,"\n") + + ## cat("est1: "); print(eff.est1 <- ida (x, y1, mcov,myCPDAG, method="local")) + ## cat("est2: "); print(eff.est2 <- ida (x, y2, mcov,myCPDAG, method="local")) + ## cat("est3: "); print(eff.est3 <- ida (x, y3, mcov,myCPDAG, method="local")) + eff.est1 <- ida (x, y1, mcov,myCPDAG, method="local") + eff.est2 <- ida (x, y2, mcov,myCPDAG, method="local") + eff.est3 <- ida (x, y3, mcov,myCPDAG, method="local") + eff.estF <- idaFast(x, c(y1,y2,y3), mcov,myCPDAG) + cat("estF:\n"); print(eff.estF) + + stopifnot(exprs = { + all.equal.1(eff.true1, eff.est1) + all.equal.1(eff.true2, eff.est2) + all.equal.1(eff.true3, eff.est3) + all.equal(unname(eff.estF), tol = 3e-15, # even tol=0 works on some + rbind(eff.est1, eff.est2, eff.est3, deparse.level=0L)) + }) + } i=1 --== estF: [,1] 9 6.465604e-01 4 0.000000e+00 2 9.627851e-17 i=2 --== estF: [,1] 4 0.000000e+00 6 0.000000e+00 8 -4.129599e-16 i=3 --== estF: [,1] [,2] [,3] 4 0.3053817 0.3053817 0.0000000 9 0.6810269 0.5676975 0.5675770 2 0.4266329 0.0000000 0.4266329 i=4 --== estF: [,1] [,2] 2 0.1670084 -2.905771e-17 10 1.5693973 8.192260e-01 1 0.4994495 0.000000e+00 i=5 --== estF: [,1] 9 0.9693677 2 0.0000000 8 0.5319225 i=6 --== estF: [,1] 9 6.395355e-17 7 0.000000e+00 1 0.000000e+00 i=7 --== estF: [,1] [,2] 3 0.1252467 0 6 0.1110212 0 4 0.2063969 0 i=8 --== estF: [,1] 4 -1.083936e-16 3 -1.342016e-16 6 -1.032320e-16 i=9 --== estF: [,1] 5 -4.587984e-17 9 -1.720494e-17 1 -9.175968e-17 i=10 --== estF: [,1] 1 0.000000 9 0.494396 4 0.000000 i=11 --== estF: [,1] [,2] [,3] [,4] [,5] 6 0.3435164 0.2833413 0.3435164 0.34351639 0.3435164 8 0.8116762 0.7946965 0.7595386 0.63218699 0.6290513 4 0.4274092 0.4274092 0.0000000 0.08282029 0.0000000 i=12 --== estF: [,1] [,2] [,3] 5 0.0000000 0 0.0000000 8 0.0000000 0 0.0000000 3 0.6774988 0 0.6774988 i=13 --== estF: [,1] 3 0.000000e+00 10 9.348925e-17 6 0.000000e+00 i=14 --== estF: [,1] 8 0.6559217 7 0.4776569 2 0.0000000 i=15 --== estF: [,1] 10 1.097852e-01 2 0.000000e+00 9 -3.516183e-17 i=16 --== estF: [,1] 10 -9.712599e-16 3 -2.869595e-16 2 1.530451e-16 i=17 --== estF: [,1] 6 0.8930411 10 2.1742263 1 0.0000000 i=18 --== estF: [,1] [,2] 6 0.5209794 4.555425e-01 2 0.4886325 0.000000e+00 1 0.1771921 3.997157e-17 i=19 --== estF: [,1] 5 0.000000e+00 8 9.281987e-01 10 -2.426762e-17 i=20 --== estF: [,1] 2 -3.623091e-17 6 0.000000e+00 1 8.453879e-17 i=21 --== estF: [,1] 7 -3.354582e-17 6 0.000000e+00 5 0.000000e+00 i=22 --== estF: [,1] 3 0.0000000 6 0.5263106 5 0.8723258 i=23 --== estF: [,1] 5 -2.293252e-17 9 4.350878e-01 3 0.000000e+00 i=24 --== estF: [,1] 5 0.4693352 4 0.6963576 3 0.0000000 i=25 --== estF: [,1] [,2] 8 0.2742938 0.2742938 4 0.9626905 0.0000000 6 0.0000000 0.0000000 i=26 --== estF: [,1] 8 -7.146787e-17 2 0.000000e+00 6 -7.146787e-17 i=27 --== estF: [,1] [,2] [,3] [,4] 1 0.1011486 0.0000000 -0.36514341 0.0000000 8 0.2830218 0.2506978 0.05199649 0.1414599 3 1.0250794 0.9588791 0.00000000 0.0000000 i=28 --== estF: [,1] [,2] [,3] [,4] 3 0.5372195 0.000000e+00 0.2824036 0.000000e+00 1 0.2327061 5.868054e-17 0.1223281 2.934027e-17 6 0.6497455 1.283729e-01 0.4024462 1.283729e-01 i=29 --== estF: [,1] [,2] 1 0.44990635 0.00000e+00 8 0.14935833 0.00000e+00 6 0.04571869 9.96617e-18 i=30 --== estF: [,1] 3 0.000000e+00 8 8.295538e-01 1 1.891401e-17 i=31 --== estF: [,1] 4 0.000000e+00 7 0.000000e+00 9 6.206343e-17 i=32 --== estF: [,1] 4 2.470293e-17 1 0.000000e+00 2 0.000000e+00 i=33 --== estF: [,1] 9 1.127628e-16 6 0.000000e+00 1 0.000000e+00 i=34 --== estF: [,1] [,2] 5 -1.259886e-16 -5.599494e-17 3 0.000000e+00 0.000000e+00 10 1.544573e+00 1.471844e+00 i=35 --== estF: [,1] 7 5.457553e-17 9 -1.571374e-16 4 1.257099e-16 i=36 --== estF: [,1] 10 0.9010311 7 0.9663906 1 0.0000000 i=37 --== estF: [,1] [,2] [,3] [,4] 8 0.8612670 0.86126699 0.8612670 8.612670e-01 5 0.0000000 0.00000000 0.0000000 0.000000e+00 10 0.4386544 0.05490201 0.3203004 5.551115e-17 i=38 --== estF: [,1] 2 -2.391233e-16 6 0.000000e+00 9 2.690137e-16 i=39 --== estF: [,1] [,2] 9 0.91109264 0.637586 5 0.06550116 0.000000 4 0.00000000 0.000000 i=40 --== estF: [,1] 1 7.409648e-17 7 9.171380e-01 9 8.729292e-01 i=41 --== estF: [,1] 1 0.0000000 5 0.8573454 7 0.6084259 i=42 --== estF: [,1] 5 3.680550e-17 6 1.472220e-16 9 9.195516e-01 i=43 --== estF: [,1] [,2] [,3] [,4] [,5] 6 0.9647508 0.7021345 0.0000000 0.9647508 0.0000000 1 0.4723034 0.0000000 0.1767421 0.4723034 0.0000000 9 0.4852022 0.3444908 0.1474874 0.4852022 0.1254846 i=44 --== estF: [,1] 3 0 8 0 2 0 i=45 --== estF: [,1] [,2] 2 0.0000000 0.0000000 8 0.1755336 0.0774081 7 0.2271521 0.1969269 i=46 --== estF: [,1] 2 0.0000000 3 0.0000000 10 0.1626532 i=47 --== estF: [,1] 2 0.0000000 7 0.7954480 3 0.6393763 i=48 --== estF: [,1] [,2] [,3] [,4] 8 0.5993879 0.3360035 0.2633844 0.5993879 9 0.6466175 0.5222904 0.1924634 0.5784813 10 2.8610582 2.3601838 1.5336939 2.8111615 i=49 --== estF: [,1] [,2] 9 1.040418 1.040418 2 0.000000 0.000000 3 0.000000 0.000000 i=50 --== estF: [,1] [,2] [,3] [,4] [,5] [,6] 4 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 9 0.5676103 0.4913266 0.3716894 0.5205935 0.3919242 0.4913266 2 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 > > proc.time() user system elapsed 2.17 0.15 2.31