tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle-model_bayes-lr.R

In npetraco/CRFutil: Handy Utility Functions for Use with CRF and gRbase packages

library(CRFutil)
library(rstan)
#library(shinystan)
#library(coda)

# Load samps generated by true model:
fpth <- "/home/npetraco/codes/R/CRFutil/tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle_data/"
load(paste0(fpth,"triangle_samp.RData"))
head(samps)

# Instantiate an empty model for the Bayesian fit:
# Fully connected
grphf <- ~1:2+1:3+2:3
adj <- ug(grphf, result="matrix")
# Check the graph:
gp <- ug(grphf, result = "graph")
dev.off()
iplot(gp)

f0 <- function(y){ as.numeric(c((y==1),(y==2)))}
n.states <- 2

bayes.lrm <- make.crf(adj, n.states)
bayes.lrm <- make.features(bayes.lrm)
bayes.lrm <- make.par(bayes.lrm, 6)
bayes.lrm$node.par[1,1,] <- 1
bayes.lrm$node.par[2,1,] <- 2
bayes.lrm$node.par[3,1,] <- 3

bayes.lrm$edge.par[[1]][1,1,1] <- 4
bayes.lrm$edge.par[[1]][2,2,1] <- 4
bayes.lrm$edge.par[[2]][1,1,1] <- 5
bayes.lrm$edge.par[[2]][2,2,1] <- 5
bayes.lrm$edge.par[[3]][1,1,1] <- 6
bayes.lrm$edge.par[[3]][2,2,1] <- 6
bayes.lrm$edges

#----------

Delta.alpha.info <- delta.alpha(crf = bayes.lrm, samples = samps, printQ = F)
Delta.alpha <- Delta.alpha.info$Delta.alpha

y <-c(samps[,1], samps[,2], samps[,3])
y[which(y==2)] <- 0
y

dat <- list(
  N=length(y),
  X1=Delta.alpha[,1],
  X2=Delta.alpha[,2],
  X3=Delta.alpha[,3],
  X4=Delta.alpha[,4],
  X5=Delta.alpha[,5],
  X6=Delta.alpha[,6],
  y = y
)


model.c <- stanc(file = "/home/npetraco/codes/R/CRFutil/tests/regression_tests/logistic_model1.stan", model_name = 'model')
rstan_options(auto_write = TRUE)
sm      <- stan_model(stanc_ret = model.c, verbose = T)

# Run the model
options(mc.cores = parallel::detectCores())
bfit <- sampling(sm, data = dat,
                 #control = list(adapt_delta = 0.85),
                 iter=20000,
                 thin = 10,
                 chains = 8)
print(bfit)

beta <- extract(bfit,"beta")[[1]]
dim(beta)


#
bayes.lrm$par <- c(median(beta[,1]), median(beta[,2]), median(beta[,3]), median(beta[,4]), median(beta[,5]), median(beta[,6]))

out.pot2 <- make.pots(parms = bayes.lrm$par,  crf = bayes.lrm,  rescaleQ = T, replaceQ = T)
bayes.lrm$node.pot
bayes.lrm$edge.pot

pot.info <- make.gRbase.potentials(bayes.lrm, node.names = gp@nodes, state.nmes = c("1","2"))
pot.info

gR.dist.info    <- distribution.from.potentials(pot.info$node.potentials, pot.info$edge.potentials)
#logZ            <- gR.dist.info$logZ
joint.dist.info <- as.data.frame(as.table(gR.dist.info$state.probs))
joint.dist.info

colnames(joint.dist.info)[4] <- c("Bayes.lr.Freq")
bayes.lr.dist.info <- joint.dist.info
save(bayes.lr.dist.info, file = paste0(fpth,"triangle_bayes-lr_dist.RData"))
bayes.lr.dist.info