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

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

# Triangle model log linear
# Bayes Poisson direct via features model matrix and Stan

# Using the output coefficients, does CRF give logZ commensurate with intercept??

#library(MASS)
library(rstan)
library(shinystan)
library(CRFutil)

# Triangle model
grphf <- ~X.1:X.2+X.1:X.3+X.2:X.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

llm2 <- make.empty.field(adj.mat = adj, parameterization.typ = "standard")

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

# Make state count freq table from samples and compute frequencies of all possible state configs.
# State configs not observed will have 0 freq
ftab <- data.frame(ftable(data.frame(samps)))
X.all <- ftab[,1:ncol(samps)]
freqs <- ftab[,ncol(ftab)]

X.all
class(X.all)
freqs

# Model Matrix with respect to graph ????
M <- compute.model.matrix(configs=X.all, edges.mat=llm2$edges, node.par = llm2$node.par, edge.par = llm2$edge.par, ff = f0)
M

# Fit log(lam) = alp + M theta
model.c <- stanc(file = paste0(fpth,"vanalla.poisson.regression2a.stan"), model_name = 'model')
sm <- stan_model(stanc_ret = model.c, verbose = T)

dat <- list(
  p = ncol(M),
  N = nrow(M),
  y = freqs,
  Mmodl = M
)
#dat

bfit <- sampling(sm, data = dat, iter=10000, thin = 4, chains = 4)
bfit
launch_shinystan(bfit)

theta <- extract(bfit, "theta")[[1]]
alpha <- extract(bfit, "alpha")[[1]]

hist(alpha)
#N <- nrow(M)
N <- nrow(samps)
logZ <- log(N) - alpha
hist(logZ)

lambda <- extract(bfit, "lambda")[[1]]
colMeans(lambda)/N      # Approx config probs
sum(colMeans(lambda)/N) # Add about to 1?



th <- colMeans(theta)
llm2$par <- th
llm2$par

out.pot <- make.pots(parms = llm2$par,  crf = llm2,  rescaleQ = F, replaceQ = T)
log(llm2$node.pot)
llm2$edge.pot

pot.info.llm2.model        <- make.gRbase.potentials(llm2, node.names = c("X.1", "X.2", "X.3"), state.nmes = c("1","2"))
gR.dist.info.llm2.model    <- distribution.from.potentials(pot.info.llm2.model$node.potentials, pot.info.llm2.model$edge.potentials)
logZ.llm2.model            <- gR.dist.info.llm2.model$logZ

logZ.llm2.model

joint.dist.info.llm2.model <- as.data.frame(as.table(gR.dist.info.llm2.model$state.probs))
joint.dist.info.llm2.model <- joint.dist.info.llm2.model[,c(2,3,1,4)]
joint.dist.info.llm2.model[,4] <- joint.dist.info.llm2.model[,4] * 100
joint.dist.info.llm2.model
#    X1 X2 X3 bayes.lrm.cp mle.lrm.cp bayes.llm.cp* bayes.llm2 bayes.llm2.lams true.model.cp
# 1  1  1  1         10.3       10.5          16.6        19.9            19.9        18.4
# 2  1  1  2          1.1        1.5          10.7         8.3             8.3         8.0
# 3  2  1  1         22.0       21.5          14.3        14.8            14.8        15.6
# 4  2  1  2         23.1       22.5          17.1        21.0            21.0        22.6
# 5  1  2  1          9.1        9.5           7.9         4.5             4.5         4.2
# 6  1  2  2         10.3       10.5          14.1        14.3            14.3        14.9
# 7  2  2  1          2.0        2.5           4.5         1.5             1.5         1.3
# 8  2  2  2         22.0       21.5          14.8        15.8            15.8        15.0

cbind(X.all, colMeans(lambda)/N*100)

config.probs.from.lam <- lambda/N*100
hist(config.probs.from.lam[,1], main=paste("X = ",X.all[1,1],X.all[1,2],X.all[1,3]), xlab="Config Prob (%)")
hist(config.probs.from.lam[,2], main=paste("X = ",X.all[2,1],X.all[2,2],X.all[2,3]), xlab="Config Prob (%)")
hist(config.probs.from.lam[,3], main=paste("X = ",X.all[3,1],X.all[3,2],X.all[3,3]), xlab="Config Prob (%)")
hist(config.probs.from.lam[,4], main=paste("X = ",X.all[4,1],X.all[4,2],X.all[4,3]), xlab="Config Prob (%)")
hist(config.probs.from.lam[,5], main=paste("X = ",X.all[5,1],X.all[5,2],X.all[5,3]), xlab="Config Prob (%)")
hist(config.probs.from.lam[,6], main=paste("X = ",X.all[6,1],X.all[6,2],X.all[6,3]), xlab="Config Prob (%)")
hist(config.probs.from.lam[,7], main=paste("X = ",X.all[7,1],X.all[7,2],X.all[7,3]), xlab="Config Prob (%)")
hist(config.probs.from.lam[,8], main=paste("X = ",X.all[8,1],X.all[8,2],X.all[8,3]), xlab="Config Prob (%)")