tests/regression_tests/hojsgaard_model_tests/triangle_model/triangle-model_mle-lr.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


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

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

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

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

#      FIX
# Ma <- glm(y ~ Delta.alpha[,1] + Delta.alpha[,2] + Delta.alpha[,3] +
#             Delta.alpha[,4] + Delta.alpha[,5] + Delta.alpha[,6] -1,
#           family=binomial(link="logit"))
Ma <- glm(y ~ Delta.alpha -1, family=binomial(link="logit"))
summary(Ma)


# c(mean(beta[,1]), mean(beta[,2]), mean(beta[,3]), mean(beta[,4]), mean(beta[,5]), mean(beta[,6]))
# c(median(beta[,1]), median(beta[,2]), median(beta[,3]), median(beta[,4]), median(beta[,5]), median(beta[,6]))
# true.model$par

mle.lrm$par <- as.numeric(coef(Ma))
mle.lrm$par

#out.pot2a <- make.pots(parms = mle.lrm$par,  crf = mle.lrm,  rescaleQ = T, replaceQ = T)
#mle.lrm$node.pot
#mle.lrm$edge.pot
pot.info2 <- make.gRbase.potentials(mle.lrm, node.names = gp@nodes, 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
mle.lr.dist.info <- joint.dist.info2[,c(2,3,1,4)]
colnames(mle.lr.dist.info)[4] <- c("MLE.lr.Freq")
mle.lr.dist.info
save(mle.lr.dist.info, file = paste0(fpth,"triangle_mle-lr_dist.RData"))