tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle-model_psl.R

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

library(CRFutil)

# Load samps generated by true model:
fpth <- "/home/npetraco/codes/R/CRFutil/tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle_data/"
setwd(fpth)
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


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

psl.m$edge.par[[1]][1,1,1] <- 4
psl.m$edge.par[[1]][2,2,1] <- 4
psl.m$edge.par[[2]][1,1,1] <- 5
psl.m$edge.par[[2]][2,2,1] <- 5
psl.m$edge.par[[3]][1,1,1] <- 6
psl.m$edge.par[[3]][2,2,1] <- 6
#----------

# Optimize theta: find minimum of negative log pseudo-likelihood:
# Auxiliary, gradient convenience function. Follows train.mrf in CRF:
gradient <- function(par, crf, ...) { crf$gradient }

#psl.m$par <- c(0,0)                  # Reset theta to start at 0
opt.info  <- stats::optim(         # optimize parameters
  par          = psl.m$par,          # theta
  fn           = neglogpseudolik2, # objective function
  gr           = gradient,         # grad of obj func
  crf          = psl.m,              # passed to fn/gr
  samples      = samps,            # passed to fn/gr
  conditional.energy.function.type = "feature", # passed to fn/gr
  ff           = f0,               # passed to fn/gr
  gradQ        = TRUE,             # passed to fn/gr
  update.crfQ  = TRUE,             # passed to fn/gr
  method       = "L-BFGS-B",
  control      = list(trace = 1, REPORT=1))

opt.info$convergence
opt.info$message
psl.m$gradient
psl.m$nll
psl.m$par

# Fit model to samples from the known model and obtain an estimate for w:
psl.m$node.pot
psl.m$par.stat
psl.m$par
#MRF theta: 0.1211475  0.8901192 -0.7510040 -0.4105264  0.6291559  1.0374155
#PSL theta: 0.1004964  0.8828685 -0.7548408 -0.4165777  0.6138656  1.0149626
psl.m$node.pot
out.pot <- make.pots(parms = psl.m$par,  crf = psl.m,  rescaleQ = T, replaceQ = T)
psl.m$node.pot


pot.info2 <- make.gRbase.potentials(psl.m, node.names = c("X.1","X.2","X.3"), state.nmes = c("1","2"))
pot.info2

gR.dist.info2    <- distribution.from.potentials(pot.info2$node.potentials, pot.info2$edge.potentials)
#logZ2            <- gR.dist.info2$logZ
joint.dist.info2 <- as.data.frame(as.table(gR.dist.info2$state.probs))
joint.dist.info2
psl.dist.info <- joint.dist.info2[,c(2,3,1,4)]
#colnames(psl.dist.info)[4] <- c("PSL.Freq")
psl.dist.info

# Compute (Besag) Pseudolikelihoods too: Pr(X) = Prod Pr(Xi|X/Xi)
tpsl.info <- pseudolikelihoods.from.energies2(state.space=psl.dist.info[,1:3], crf = psl.m, cond.en.form="feature", ff=f0)
tpsl.info$condtional.energies
tpsl.info$complement.condtional.energies
tpsl.info$conditional.Zs
tpsl.info$conditional.Prs
tpsl.info$complement.conditional.Prs
tpsl.info$pseudo.likelihoods

psl.dist.info <- cbind(psl.dist.info, tpsl.info$pseudo.likelihoods)
colnames(psl.dist.info)[c(4,5)] <- c("PSL.MRF.Freq","PSL.Freq")
psl.dist.info

save(psl.dist.info, file = paste0(fpth,"triangle_psl_dist.RData"))