########## R script: BobTheDAGwithData ########## # For conducting examples involving Bob the DAG # with the help of the "rjags" package in R. # Last changed: 30 AUG 2023 # Load in the rjags package and the writeModel() function: library(rjags) source("writeModel.r") # Set Monte Carlo sample sizes: nBurnin <- 10000 ; nIter <- 10000 ; nThin <- 1 # Specify P(x2=1) value: probx2eq1 <- 0.61 # Set the data: x1Value <- 0 # Set up input lists for rjags: allData <- list(probx2eq1=probx2eq1,x1=x1Value) initValues <- list(list(x2=0)) # Specify model in JAGS: BobTheDAGModel <- function() { x2 ~ dbern(probx2eq1) probx1eq1CondOnx2 <- (1*x2 + 2)/(3*x2 + 7) x1 ~ dbern(probx1eq1CondOnx2) } modelFileName <- file.path(getwd(),"BobTheDAGModel.txt") writeModel(BobTheDAGModel,modelFileName) # Generate the Monte Carlo samples: currModel<-jags.model("BobTheDAGModel.txt",data=allData) update(currModel,n.iter=nBurnin) currMCMC <- jags.samples(currModel, variable.names="x2", n.iter=nIter,thin=nThin) x2MC <- as.numeric(currMCMC$x2) # Obtain Monte Carlo approximations to probability functions: px2MC <- table(x2MC)/length(x2MC) # Obtain exact values: numerCoefx2 <- 1 ; numerConst <- 2 denomCoefx2 <- 3 ; denomConst <- 7 wp <- probx2eq1 bN <- numerCoefx2 ; aN <- numerConst bD <- denomCoefx2 ; aD <- denomConst if (x1Value==0) { numer <- wp*(1 - ((aN + bN )/(aD + bD))) denom <- (1-wp)*((aD-aN)/aD) + wp*(((aD-aN)+(bD-bN))/(aD+bD)) px2eq1givx1eq0 <- numer/denom px2givenx1 <- c(1-px2eq1givx1eq0,px2eq1givx1eq0) } if (x1Value==1) { numer <- wp*((aN + bN )/(aD + bD)) denom <- (1-wp)*(aN/aD) + wp*((aN+bN)/(aD+bD)) px2eq1givx1eq1 <- numer/denom px2givenx1 <- c(1-px2eq1givx1eq1,px2eq1givx1eq1) } outMatx2 <- rbind(px2givenx1,px2MC) outMatx2 <- as.data.frame(outMatx2) dimnames(outMatx2)[[1]] <- c("exact","Monte Carlo") dimnames(outMatx2)[[2]] <- c("x2=0","x2=1") cat("\n\n\n") cat("Results for p(x2|x1) when value of x1 is ",x1Value,":\n",sep="") print(round(outMatx2,4)) cat("\n\n\n") # Plot exact and Monte Carlo probability functions: colEX <- "salmon" ; colMC <- "blue" par(mfrow=c(2,2),mai=c(0.4,0.7,0.5,0.08)) plot(0,0,type="n",xlim=c(0,1),ylim=c(0,1),xaxt="n",yaxt="n", xlab="",ylab="",bty="o") text(0.5,0.8,"Bob the",cex=2,col="darkgreen") text(0.5,0.5,"DAG",cex=2,col="darkgreen") text(0.5,0.2,"with data",cex=2,col="darkgreen") plot(0,0,type="n",xlim=c(0,1),ylim=c(0,1),xaxt="n",yaxt="n", xlab="",ylab="",bty="o") text(0.5,0.7,"DATA:",col="red",cex=2) if (x1Value==0) textString <- expression(paste(x[1],"=0")) if (x1Value==1) textString <- expression(paste(x[1],"=1")) text(0.5,0.4,textString,col="red",cex=2) plot(0,0,type="n",xlim=c(0,1),ylim=c(0,1),xaxt="n",yaxt="n", xlab="",ylab="",bty="o") legend(0.05,0.6,legend=c("exact","Monte Carlo"), col=c(colEX,colMC),lwd=c(2,7),cex=1,bty="n") plot(px2MC,bty="n",lwd=7,col=colMC,ylab="probability", ylim=c(0,max(c(px2MC,px2givenx1))),xlab="", main=expression(paste("p(",x[2],"|",x[1],")"))) lines(rep(0,2),c(0,px2givenx1[1]),col=colEX,lwd=2) lines(rep(1,2),c(0,px2givenx1[2]),col=colEX,lwd=2) abline(0,0,col="navy") ########## End of BobTheDAGwithData ##########