########## R script: NathanielTheDAGwithData ########## # For conducting examples involving Nathaniel the DAG # with data, with the help of the "rjags" package in R. # Last changed: 09 AUG 2022 # Load in the rjags package and the writeModel() function: library(rjags) source("writeModel.r") nBurnin <- 10000 ; nIter <- 10000 ; nThin <- 100 pX6 <- c(0.11,0.35,0.07,0.41,0.06) pX7 <- c(0.21,0.28,0.10,0.18,0.23) pX8 <- c(0.34,0.04,0.25,0.17,0.20) # Set data: x1Value <- 4 ; x2Value <- 2 # Do Monte Carlo approximations: allData <- list(pX6=pX6,pX7=pX7,pX8=pX8,X1=x1Value,X2=x2Value) initValues <- list(list(X3=3,X4=3,X5=3,X6=3,X7=3,X8=3)) # Specify model in JAGS: NathanielTheDAGModel <- function() { X6 ~ dcat(pX6[]) X7 ~ dcat(pX7[]) X8 ~ dcat(pX8[]) # Determine p(x4): for (x4 in 1:5) { psi46[x4] <- abs(x4 - X6) + 0.134 } for (x4 in 1:5) { pX4[x4] <- psi46[x4]/sum(psi46[]) } X4 ~ dcat(pX4[]) # Determine p(x5): for (x5 in 1:5) { psi578[x5] <- sqrt((x5 + 0.7*X7)/(0.9 + X8)) } for (x5 in 1:5) { pX5[x5] <- psi578[x5]/sum(psi578[]) } X5 ~ dcat(pX5[]) # Determine p(x2): for (x2 in 1:5) { psi245[x2] <- pow(x2 + 0.7*X4 + X5,0.2) } for (x2 in 1:5) { pX2[x2] <- psi245[x2]/sum(psi245[]) } X2 ~ dcat(pX2[]) # Determine p(x3): for (x3 in 1:5) { psi35[x3] <- log(cosh(x3 + X5 + 2)) } for (x3 in 1:5) { pX3[x3] <- psi35[x3]/sum(psi35[]) } X3 ~ dcat(pX3[]) # Determine p(x1): for (x1 in 1:5) { psi13[x1] <- pow(sinh(x1 + 0.3*X3 + 1.8),2) } for (x1 in 1:5) { pX1[x1] <- psi13[x1]/sum(psi13[]) } X1 ~ dcat(pX1[]) } modelFileName <- file.path(getwd(),"NathanielTheDAGModel.txt") writeModel(NathanielTheDAGModel,modelFileName) # Run rjags: currModel<-jags.model("NathanielTheDAGModel.txt",data=allData) update(currModel,n.iter=nBurnin) currMCMC <- jags.samples(currModel, variable.names=c("X3","X4","X5","X6","X7","X8"), n.iter=nIter,thin=nThin) # Extract MC samples: X3MC <- as.numeric(currMCMC$X3) X4MC <- as.numeric(currMCMC$X4) X5MC <- as.numeric(currMCMC$X5) X6MC <- as.numeric(currMCMC$X6) X7MC <- as.numeric(currMCMC$X7) X8MC <- as.numeric(currMCMC$X8) # Obtain MC approximations to probability functions: p3MC <- table(X3MC)/length(X3MC) p4MC <- table(X4MC)/length(X4MC) p5MC <- table(X5MC)/length(X5MC) p6MC <- table(X6MC)/length(X6MC) p7MC <- table(X7MC)/length(X7MC) p8MC <- table(X8MC)/length(X8MC) # Do exact calculations, starting with setting of # all potential functions: psi13 <- function(x1,x3) return(sinh(x1 + 0.3*x3 + 1.8)^2) psi245 <- function(x2,x4,x5) return((x2 + 0.7*x4 + x5)^0.2) psi35 <- function(x3,x5) return(log(cosh(x3 + x5+2))) psi46 <- function(x4,x6) return(abs(x4 - x6) + 0.134) psi578 <- function(x5,x7,x8) return(sqrt((x5 + 0.7*x7)/(0.9 + x8))) # Obtain p(x1|x3): pX1condX3 <- function(x1,x3) { normFac <- 0 for (x1Dash in 1:5) normFac <- normFac + psi13(x1Dash,x3) return(psi13(x1,x3)/normFac) } # Obtain p(x2|x4,x5): pX2condX4X5 <- function(x2,x4,x5) { normFac <- 0 for (x2Dash in 1:5) normFac <- normFac + psi245(x2Dash,x4,x5) return(psi245(x2,x4,x5)/normFac) } # Obtain p(x3|x5): pX3condX5 <- function(x3,x5) { normFac <- 0 for (x3Dash in 1:5) normFac <- normFac + psi35(x3Dash,x5) return(psi35(x3,x5)/normFac) } # Obtain p(x4|x6): pX4condX6 <- function(x4,x6) { normFac <- 0 for (x4Dash in 1:5) normFac <- normFac + psi46(x4Dash,x6) return(psi46(x4,x6)/normFac) } # Obtain p(x5|x7,x8): pX5condX7X8 <- function(x5,x7,x8) { normFac <- 0 for (x5Dash in 1:5) normFac <- normFac + psi578(x5Dash,x7,x8) return(psi578(x5,x7,x8)/normFac) } # Obtain p(x3|x1Value,x2Value): pX3condX1X2UN <- rep(0,5) for (x3 in 1:5) for (x4 in 1:5) for (x5 in 1:5) for (x6 in 1:5) for (x7 in 1:5) for (x8 in 1:5) pX3condX1X2UN[x3] <- pX3condX1X2UN[x3] + (pX1condX3(x1Value,x3) *pX2condX4X5(x2Value,x4,x5)*pX3condX5(x3,x5)*pX4condX6(x4,x6) *pX5condX7X8(x5,x7,x8)*pX6[x6]*pX7[x7]*pX8[x8]) pX3condX1X2 <- pX3condX1X2UN/sum(pX3condX1X2UN) # Obtain p(x4|x1Value,x2Value): pX4condX1X2UN <- rep(0,5) for (x4 in 1:5) for (x3 in 1:5) for (x5 in 1:5) for (x6 in 1:5) for (x7 in 1:5) for (x8 in 1:5) pX4condX1X2UN[x4] <- pX4condX1X2UN[x4] + (pX1condX3(x1Value,x3) *pX2condX4X5(x2Value,x4,x5)*pX3condX5(x3,x5)*pX4condX6(x4,x6) *pX5condX7X8(x5,x7,x8)*pX6[x6]*pX7[x7]*pX8[x8]) pX4condX1X2 <- pX4condX1X2UN/sum(pX4condX1X2UN) # Obtain p(x5|x1Value,x2Value): pX5condX1X2UN <- rep(0,5) for (x5 in 1:5) for (x3 in 1:5) for (x4 in 1:5) for (x6 in 1:5) for (x7 in 1:5) for (x8 in 1:5) pX5condX1X2UN[x5] <- pX5condX1X2UN[x5] + (pX1condX3(x1Value,x3) *pX2condX4X5(x2Value,x4,x5)*pX3condX5(x3,x5)*pX4condX6(x4,x6) *pX5condX7X8(x5,x7,x8)*pX6[x6]*pX7[x7]*pX8[x8]) pX5condX1X2 <- pX5condX1X2UN/sum(pX5condX1X2UN) # Obtain p(x6|x1Value,x2Value): pX6condX1X2UN <- rep(0,5) for (x6 in 1:5) for (x3 in 1:5) for (x4 in 1:5) for (x5 in 1:5) for (x7 in 1:5) for (x8 in 1:5) pX6condX1X2UN[x6] <- pX6condX1X2UN[x6] + (pX1condX3(x1Value,x3) *pX2condX4X5(x2Value,x4,x5)*pX3condX5(x3,x5)*pX4condX6(x4,x6) *pX5condX7X8(x5,x7,x8)*pX6[x6]*pX7[x7]*pX8[x8]) pX6condX1X2 <- pX6condX1X2UN/sum(pX6condX1X2UN) # Obtain p(x7|x1Value,x2Value): pX7condX1X2UN <- rep(0,5) for (x7 in 1:5) for (x3 in 1:5) for (x4 in 1:5) for (x5 in 1:5) for (x6 in 1:5) for (x8 in 1:5) pX7condX1X2UN[x7] <- pX7condX1X2UN[x7] + (pX1condX3(x1Value,x3) *pX2condX4X5(x2Value,x4,x5)*pX3condX5(x3,x5)*pX4condX6(x4,x6) *pX5condX7X8(x5,x7,x8)*pX6[x6]*pX7[x7]*pX8[x8]) pX7condX1X2 <- pX7condX1X2UN/sum(pX7condX1X2UN) # Obtain p(x8|x1Value,x2Value): pX8condX1X2UN <- rep(0,5) for (x8 in 1:5) for (x3 in 1:5) for (x4 in 1:5) for (x5 in 1:5) for (x6 in 1:5) for (x7 in 1:5) pX8condX1X2UN[x8] <- pX8condX1X2UN[x8] + (pX1condX3(x1Value,x3) *pX2condX4X5(x2Value,x4,x5)*pX3condX5(x3,x5)*pX4condX6(x4,x6) *pX5condX7X8(x5,x7,x8)*pX6[x6]*pX7[x7]*pX8[x8]) pX8condX1X2 <- pX8condX1X2UN/sum(pX8condX1X2UN) # Do plots to compare Monte Carlo and exact results: colEX <- "salmon" ; colMC <- "blue" par(mfrow=c(3,3),mai=c(0.5,0.8,0.05,0.1)) 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,"Nathaniel",cex=2,col="darkgreen") text(0.5,0.5,"the 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.8,"DATA:",col="red",cex=2) if (x1Value==1) textString1 <- expression(paste(x[1],"=1")) if (x1Value==2) textString1 <- expression(paste(x[1],"=2")) if (x1Value==3) textString1 <- expression(paste(x[1],"=3")) if (x1Value==4) textString1 <- expression(paste(x[1],"=4")) if (x1Value==5) textString1 <- expression(paste(x[1],"=5")) if (x2Value==1) textString2 <- expression(paste(x[2],"=1")) if (x2Value==2) textString2 <- expression(paste(x[2],"=2")) if (x2Value==3) textString2 <- expression(paste(x[2],"=3")) if (x2Value==4) textString2 <- expression(paste(x[2],"=4")) if (x2Value==5) textString2 <- expression(paste(x[2],"=5")) text(0.5,0.5,textString1,col="red",cex=2) text(0.5,0.3,textString2,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(p3MC,bty="n",lwd=7,col=colMC,ylab="probability", ylim=c(0,max(c(p3MC,pX3condX1X2))),xlab="", main=expression(paste("p(",x[3],"|",x[1],",",x[2],")"))) lines(as.table(pX3condX1X2),col=colEX,lwd=2) abline(0,0,col="navy") plot(p4MC,bty="n",lwd=7,col=colMC,ylab="probability", ylim=c(0,max(c(p4MC,pX4condX1X2))),xlab="", main=expression(paste("p(",x[4],"|",x[1],",",x[2],")"))) lines(as.table(pX4condX1X2),col=colEX,lwd=2) abline(0,0,col="navy") plot(p5MC,bty="n",lwd=7,col=colMC,ylab="probability", ylim=c(0,max(c(p5MC,pX5condX1X2))),xlab="", main=expression(paste("p(",x[5],"|",x[1],",",x[2],")"))) lines(as.table(pX5condX1X2),col=colEX,lwd=2) abline(0,0,col="navy") plot(p6MC,bty="n",lwd=7,col=colMC,ylab="probability", ylim=c(0,max(c(p6MC,pX6))),xlab="", main=expression(paste("p(",x[6],"|",x[1],",",x[2],")"))) lines(as.table(pX6condX1X2),col=colEX,lwd=2) abline(0,0,col="navy") plot(p7MC,bty="n",lwd=7,col=colMC,ylab="probability", ylim=c(0,max(c(p7MC,pX7))),xlab="", main=expression(paste("p(",x[7],"|",x[1],",",x[2],")"))) lines(as.table(pX7condX1X2),col=colEX,lwd=2) abline(0,0,col="navy") plot(p8MC,bty="n",lwd=7,col=colMC,ylab="probability", ylim=c(0,max(c(p8MC,pX8))),xlab="", main=expression(paste("p(",x[8],"|",x[1],",",x[2],")"))) lines(as.table(pX8condX1X2),col=colEX,lwd=2) abline(0,0,col="navy") ########## End of NathanielTheDAGwithData ##########