Commit 622d0b4b authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

paralell tempering

parent 12cf1e54
Pipeline #2779 passed with stage
in 6 seconds
export(mcmcabn,
CoupledHeatedmcmcabn,
query,
plot.mcmcabn,
summary.mcmcabn,
print.summary.mcmcabn
print.summary.mcmcabn,
print.mcmcabn
)
importFrom(gRbase, topoSortMAT)
......@@ -19,6 +21,7 @@ print.summary.mcmcabn
#methods for class mcmcabn
S3method(plot, mcmcabn)
S3method(summary, mcmcabn)
S3method(print, mcmcabn)
#method for class summary.mcmcabn
S3method(print, summary.mcmcabn)
......@@ -9,4 +9,5 @@
## mcmcabn 0.3:
* update mcmcabn to make it compatible with constraints imported from the cache of scores. Heating parameter to increase or decrease acceptance probability.
* new function CoupledHeatedmcmcabn() implementing parallel tempering
* new article
......@@ -198,14 +198,11 @@ if(FALSE){
s.proposed <- score.dag(dag.MBR,score.cache,sc)
s.current <- score.dag(dag.tmp,score.cache,sc)
A <- min(exp(( s.proposed - s.current) * (n.edges/n.edges.tilde) ), 1)
A <- min(exp( (s.proposed - s.current) * (n.edges/n.edges.tilde)), 1)
if(is.nan(A)){A <- 0}
# if((score.A)<0){score.A <- 0}
# A <- min(1, score.A)
if (runif(1)<(A**heating)) {
if (runif(1)<exponent(A, heating)) {
rejection <- 0
dag.tmp <- dag.MBR
......
......@@ -117,7 +117,7 @@ REV <- function(n.var, dag.tmp, retain, ban, max.parents, sc, score.cache, score
#A <- min(1, score.A)
#if (rbinom(n = 1, size = 1, prob = A) == 1) {
if (runif(1)<(A**heating)) {
if (runif(1)<(exponent(A,heating))) {
rejection <- 0
dag.tmp <- dag.M.tilde
score <- score.dag(dag.M.tilde,score.cache,sc)
......
......@@ -144,13 +144,13 @@ mc3 <- function(n.var, dag.tmp, retain, ban, max.parents, sc, score.cache, score
score.Gprime.scaled <- score.dag(dag.gprime,score.cache,sc)
score.G.scaled <- score.dag(dag.tmp,score.cache,sc)
alpha <- min(exp(( score.Gprime.scaled - score.G.scaled) * (n.G/n.Gprime) * (prior.Gprime/prior.G)), 1)
alpha <- min(exp( (score.Gprime.scaled - score.G.scaled) * (n.G/n.Gprime) * (prior.Gprime/prior.G)), 1)
if(!is.numeric(alpha) | is.nan(alpha)) alpha <- 0
score <- score.G
if (!is.null(dag.gprime) && runif(1)<(alpha**heating)) {
if (!is.null(dag.gprime) && runif(1)<exponent(alpha,heating)) {
dag.tmp <- dag.gprime
score <- score.Gprime
rejection <- 0
......
##############################################################################
## mcmcabn.R --- Author : Gilles Kratzer Document created: 28/11/2019 -
##-------------------------------------------------------------------------
## Coupled heated MCMC procedure
##-------------------------------------------------------------------------
CoupledHeatedmcmcabn <- function(score.cache = NULL, score = "mlik", data.dists = NULL, max.parents = 1, mcmc.scheme = c(100, 1000,
1000), seed = 42, verbose = FALSE, start.dag = NULL, prior.dag = NULL, prior.lambda = NULL, prob.rev = 0.05, prob.mbr = 0.05, heating = 1, n.chains = 4,
prior.choice = 2) {
#################################################### Tests
.tests.mcmcabn(score.cache, data.dists, max.parents, mcmc.scheme, seed, verbose, start.dag, prior.dag, prior.lambda,
prob.rev, prob.mbr, prior.choice,heating)
if(n.chains<2){stop("Coupled heated MCMC runs should have two chains at minimum.")}
## end of tests
## format
if (is.character(score))
score <- tolower(score)
score <- c("bic", "aic", "mdl", "mlik")[pmatch(score, c("bic", "aic", "mdl", "mlik"))]
if (score %!in% c("bic", "aic", "mdl", "mlik"))
stop("A method should be provided and be one of: bic, aic, mdl or mlik.")
if (is.matrix(prior.dag)) {
prior.choice <- 3
}
if (is.matrix(prior.dag) && is.null(prior.lambda)) {
prior.lambda <- 1
}
n.var <- length(data.dists)
prob.mc3 <- 1 - (prob.rev + prob.mbr)
## output
out.dags <- array(data = NA, dim = c(n.var, n.var, mcmc.scheme[1] + 1))
out.scores <- array(data = NA, dim = c(mcmc.scheme[1] + 1))
out.scores.coupled <- array(data = NA, dim = c(mcmc.scheme[1] + 1, n.chains))
out.alpha <- array(data = NA, dim = c(mcmc.scheme[1] + 1))
out.method <- array(data = NA, dim = c(mcmc.scheme[1] + 1))
out.rejection <- array(data = NA, dim = c(mcmc.scheme[1] + 1))
heating.para <- array(data = NA, dim = c(n.chains))
## seeding
set.seed(seed = seed)
## structural restriction
retain <- score.cache$dag.retained
ban <- score.cache$dag.banned
if(is.null(retain)){retain <- matrix(data = 0, nrow = n.var, ncol = n.var)}
if(is.null(ban)){ban <- matrix(data = 0, nrow = n.var, ncol = n.var)}
## Initializing matrix
dag.tmp <- matrix(data = 0, nrow = n.var, ncol = n.var)
## start zero matrix if(!is.null(start.dag)){ colnames(dag.tmp) <- rownames(dag.tmp)<-name.dag<-sample(1:n.var) }
## start random matrix
if (is.null(start.dag)) {
start.dag <- "random"
}
if (is.character(start.dag) && start.dag == "random") {
vec.tmp <- c(rep(1, max.parents), rep(0, 2 * n.var - max.parents))
for (lines in 1:(n.var - 1)) {
dag.tmp[1 + lines, 1:lines] <- sample(vec.tmp)[1:lines]
}
colnames(dag.tmp) <- rownames(dag.tmp) <- name.dag <- sample(1:n.var)
dag.tmp <- dag.tmp[order(name.dag), order(name.dag)]
#check constraints
dag.tmp <- dag.tmp + retain
dag.tmp <- dag.tmp*(1-ban)
## run a series of checks on the DAG passed
dag.tmp <- abs(dag.tmp)
## consistency checks
diag(dag.tmp) <- 0
dag.tmp[dag.tmp > 0] <- 1
}
if (is.character(start.dag) && start.dag == "hc"){
start.dag <- searchHeuristic(score.cache = score.cache,
score = score,
num.searches = 100,
max.steps = 500,
seed = seed,
verbose = verbose,
start.dag = NULL,
dag.retained = NULL,
dag.banned = NULL,
algo = "hc",
tabu.memory = 10,
temperature = 0.9)
start.dag <- start.dag[["dags"]][[which.max(x = unlist(start.dag$scores))]]
}
if (is.matrix(start.dag)) {
start.dag[start.dag != 0] <- 1
diag(start.dag) <- 0
dag.tmp <- start.dag
colnames(dag.tmp) <- rownames(dag.tmp) <- 1:n.var
}
## scores
if (score %in% c("bic", "aic", "mdl")) {
tmp <- score.cache$node.defn[, as.numeric(colnames(dag.tmp))]
colnames(tmp) <- as.numeric(colnames(dag.tmp))
sc <- cbind(tmp, -score.cache[[score]])
}
if (score == "mlik") {
tmp <- score.cache$node.defn[, as.numeric(colnames(dag.tmp))]
colnames(tmp) <- as.numeric(colnames(dag.tmp))
sc <- cbind(tmp, score.cache[[score]])
}
## scoring init
score.init <- score.dag(dag.tmp,score.cache,sc)
## out
out.dags[, , 1] <- dag.tmp
out.scores[1] <- score <- score.init
out.scores.coupled[1,] <- rep(score.init, n.chains)
out.alpha[1] <- alpha <- 0
out.method[1] <- "MC3"
out.rejection[1] <- 0
for (j in 1:n.chains) {
heating.para[j] <- 1/(1+heating*(j-1))
}
## start mcmc search
if (verbose)
cat("Start MCMC burn in \n")
j <- 1
tmp.dag <- array(data = NA, dim = c(n.var, n.var, n.chains))
tmp.scores <- array(data = NA, dim = c(n.chains))
tmp.alpha <- array(data = NA, dim = c(n.chains))
tmp.rejection <- array(data = NA, dim = c(n.chains))
tmp.method <- array(data = NA, dim = c(n.chains))
#init
tmp.dag[,,] <- dag.tmp
tmp.scores[] <- score.init
while (j <= mcmc.scheme[3]) {
## method choice:
for (c in 1:n.chains) {
heating <- heating.para[c]
dag.tmp <- tmp.dag[,,c]
score <- tmp.scores[c]
#print(score)
method.choice <- sample(x = c("MC3", "REV", "MBR"), size = 1, prob = c(prob.mc3, prob.rev, prob.mbr))
switch(method.choice, MC3 = {
out <- mc3(n.var, (dag.tmp), retain, ban, max.parents, sc, score.cache, score, prior.choice, prior.lambda, prior.dag,
verbose, heating)
tmp.dag[,,c] <- out$dag.tmp
tmp.scores[c] <- out$score
}, REV = {
if (verbose) {
print("REV move")
}
out <- REV(n.var, (dag.tmp),retain, ban, max.parents, sc, score.cache, score, verbose, heating)
tmp.dag[,,c] <- out$dag.tmp
tmp.scores[c] <- out$score
}, MBR = {
if (verbose) {
print("MBR move")
}
out <- MBR(n.var, dag.tmp,retain, ban, max.parents, sc, score.cache, score, verbose, heating)
tmp.dag[,,c] <- out$dag.tmp
tmp.scores[c] <- out$score
})
}#eof: multiple heated chains
## swap between two randomly chosen chains
rand.chains <- sample(x = 1:n.chains,2)
alpha <- min(exponent(exp(tmp.scores[rand.chains[1]] - tmp.scores[rand.chains[2]]),heating.para[rand.chains[1]]) * exponent(exp(tmp.scores[rand.chains[2]] - tmp.scores[rand.chains[1]]),heating.para[rand.chains[2]]),1)
if(is.nan(alpha)){alpha <- 0}
if (runif(1)<(alpha)) {
tmp.dag[,,rand.chains[1]] <- tmp.dag[,,rand.chains[2]]
tmp.scores[rand.chains[1]] <- tmp.scores[rand.chains[2]]
}
j <- j + 1
} #EOF Burn in
out.scores[1] <- tmp.scores[1]
out.scores.coupled[1,] <- tmp.scores
out.dags[,,1] <- tmp.dag[,,1]
### EOF: Burn-in phase!
if (verbose)
cat("Start MCMC search \n")
for (mcmc.search in 1:mcmc.scheme[1]) {
if (verbose)
cat("processing search ...", mcmc.search, "\n")
## start blind search
j <- 0
while (j <= mcmc.scheme[2]) {
method.choice <- sample(x = c("MC3", "REV", "MBR"), size = 1, prob = c(prob.mc3, prob.rev, prob.mbr))
tmp.method[] <- method.choice
for (c in 1:n.chains) {
heating <- heating.para[c]
dag.tmp <- tmp.dag[,,c]
score <- tmp.scores[c]
switch(method.choice, MC3 = {
out <- mc3(n.var, (dag.tmp), retain, ban, max.parents, sc, score.cache, score, prior.choice, prior.lambda, prior.dag,
verbose, heating)
tmp.dag[,,c] <- out$dag.tmp
tmp.scores[c] <- out$score
#tmp.alpha[c] <- out$alpha
tmp.rejection[c] <- out$rejection
}, REV = {
if (verbose) {
print("REV move")
}
out <- REV(n.var, (dag.tmp),retain, ban, max.parents, sc, score.cache, score, verbose, heating)
tmp.dag[,,c] <- out$dag.tmp
tmp.scores[c] <- out$score
#tmp.alpha[c] <- out$alpha
tmp.rejection[c] <- out$rejection
}, MBR = {
if (verbose) {
print("MBR move")
}
out <- MBR(n.var, (dag.tmp), retain, ban, max.parents, sc, score.cache, score, verbose, heating)
tmp.dag[,,c] <- out$dag.tmp
tmp.scores[c] <- out$score
#tmp.alpha[c] <- out$alpha
tmp.rejection[c] <- out$rejection
})
}#eof: multiple heated chains
## swap between two randomly chosen chains
rand.chains <- sample(x = 1:n.chains,2)
alpha <- min(exponent(exp(tmp.scores[rand.chains[1]] - tmp.scores[rand.chains[2]]),heating.para[rand.chains[1]]) * exponent(exp(tmp.scores[rand.chains[2]] - tmp.scores[rand.chains[1]]),heating.para[rand.chains[2]]),1)
if(is.nan(alpha)){alpha <- 0}
if (runif(1)<(alpha)) {
storage.tmp1 <- tmp.dag[,,rand.chains[1]]
storage.tmp2 <-tmp.scores[rand.chains[1]]
storage.tmp3 <- tmp.method[rand.chains[1]]
tmp.dag[,,rand.chains[1]] <- tmp.dag[,,rand.chains[2]]
tmp.scores[rand.chains[1]] <- tmp.scores[rand.chains[2]]
tmp.method[rand.chains[1]] <- tmp.method[rand.chains[2]]
tmp.dag[,,rand.chains[2]] <- storage.tmp1
tmp.scores[rand.chains[2]] <- storage.tmp2
tmp.method[rand.chains[2]] <- storage.tmp3
#rejection <- 0
}
j <- j + 1
} #EOF mcmc search
out.dags[, , mcmc.search + 1] <- tmp.dag[,,1]
out.scores[mcmc.search + 1] <- tmp.scores[1]
out.alpha[mcmc.search + 1] <- alpha
out.method[mcmc.search + 1] <- tmp.method[1]
out.rejection[mcmc.search + 1] <- tmp.rejection[1]
out.scores.coupled[mcmc.search + 1,] <- tmp.scores
out.scores.coupled <- data.frame(out.scores.coupled)
names(out.scores.coupled) <- 1:n.chains
} #EOF mcmc search
out <- list(dags = out.dags, scores = out.scores, alpha = out.alpha, method = out.method, rejection = out.rejection,
iterations = mcmc.scheme[1] * (mcmc.scheme[2] + 1), thinning = mcmc.scheme[2], burnin = mcmc.scheme[3], data.dist = data.dists, heating = heating.para, score.coupled = out.scores.coupled)
class(out) <- "mcmcabn"
return(out)
} #eof
......@@ -245,3 +245,9 @@ score.dag <- function(dag,bsc.score,sc){
}
return(score)
}
##-------------------------------------------------------------------------
## Raising negative numbers to a fractional exponent
##-------------------------------------------------------------------------
exponent <- function(a, pow) (abs(a)^pow)*sign(a)
% mcmcabn.Rd ---
% Author : Gilles Kratzer
% Created on : 28.11.2019
% Last modification :
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\name{CoupledHeatedmcmcabn}
\alias{CoupledHeatedmcmcabn}
\title{Coupled Heated Structural MCMC sampler for DAGs}
\usage{
CoupledHeatedmcmcabn(score.cache = NULL,
score = "mlik",
data.dists = NULL,
max.parents = 1,
mcmc.scheme = c(100,1000,1000),
seed = 42,
verbose = FALSE,
start.dag = NULL,
prior.dag = NULL,
prior.lambda = NULL,
prob.rev = 0.05,
prob.mbr = 0.05,
heating = 1,
n.chains = 4,
prior.choice = 2)
}
\arguments{
\item{score.cache}{output from \link[abn:build_score_cache]{buildscorecache} from the \code{abn} R package.}
\item{score}{character giving which network score should be used to sample the DAGs landscape.}
\item{data.dists}{a named list giving the distribution for each node in the network, see details.}
\item{max.parents}{a constant giving the maximum number of parents allowed.}
\item{mcmc.scheme}{a sampling scheme. It is vector giving in that order: the number of returned DAGS, the number of thinned steps and length of the burn-in phase.}
\item{seed}{a non-negative integer which sets the seed.}
\item{verbose}{extra output, see output for details.}
\item{start.dag}{a DAG given as a matrix, see details for format, which can be used to provide a starting point for the structural search explicitly. Alternatively, character "random" will select a random DAG as a starting point. Character "hc" will call a hill-climber to select a DAG as a starting point.}
\item{prior.dag}{user defined prior. It should be given as a matrix where entries range from zero to one. 0.5 is non-informative for the given arc.}
\item{prior.lambda}{hyper parameter representing the strength of belief in the user-defined prior.}
\item{prob.rev}{probability of selecting a new edge reversal.}
\item{prob.mbr}{probability of selecting a Markov blanket resampling move.}
\item{heating}{a real positive number that heats up the chains. The default is one. See details}
\item{n.chains}{Number of chains to be run, see details.}
\item{prior.choice}{an integer, 1 or 2, where 1 is a uniform structural prior and 2 uses a weighted prior, see details.}
}
\description{
This function is a coupled heated structural Monte Carlo Markov Chain sampler that is equipped with two large scale MCMC moves and parallel tempering that are purposed to synergitically accelerate chain mixing.
}
\details{The procedure runs a coupled heated structural Monte Carlo Markov Chain to find the most probable posterior network (DAG). The default algorithm is based on three MCMC moves in a parallel tempering scheme: edge addition, edge deletion, and edge reversal. This algorithm is known as the (MC)^3. It is known to mix slowly and getting stuck in low probability regions. Indeed, changing of Markov equivalence region often requires multiple MCMC moves. Then large scale MCMC moves are implemented. The user can set the relative frequencies. The new edge reversal move (REV) from Grzegorczyk and Husmeier (2008) and the Markov blanket resampling (MBR) from Su and Borsuk (2016). The classical reversal move depends on the global configuration of the parents and children and fails to propose MCMC jumps that produce valid but very different DAGs in a unique move. The REV move sample globally a new set of parents. The MBR workaround applies the same idea but to the entire Markov blanket of a randomly chosen node.
The classical (MC)^3 is unbiased but inefficient in mixing. The two radical MCMC alternative moves are known to accelerate mixing without introducing biases. Those MCMC moves are computationally expensive. Then low frequencies are advised. The REV move is not necessarily ergotic, then it should not be used alone.
The parallel tempering scheme has been proposed by (Geyer, 1991). The idea is to run multiple MCMC chains, at different temperature. This will flatten the posterior probability distribution and then making jumps across low probability regions more probable. At each iteration, a swap between two randomly chosen chains is evaluated through a Metropolis like probability. THe temperature is sequentially increased in the chains folowwing
The parameter \code{start.dag} can be: \code{"random"}, \code{"hc"} or user defined. If user select \code{"random"} then a random valid DAG is selected. The routine used favourise low density structure. If \code{"hc"} (for Hill-climber: \link[abn:search_heuristic]{searchHeuristic} then a DAG is selected using 100 different searches with 500 optimization steps. A user defined DAG can be provided. It should be a named square matrix containing only zeros and ones. The DAG should be valid (i.e. acyclic).
The parameter \code{prior.choice} determines the prior used within each node for a given choice of parent combination. In Koivisto and Sood (2004) p.554, a form of prior is used, which assumes that the prior probability for parent combinations comprising of the same number of parents are all equal. Specifically, that the prior probability for parent set G with cardinality |G| is proportional to 1/[n-1 choose |G|] where there are n total nodes. Note that this favors parent combinations with either very low or very high cardinality, which may not be appropriate. This prior is used when \code{prior.choice=2}. When \code{prior.choice=1} an uninformative prior is used where parent combinations of all cardinalities are equally likely. When \code{prior.choice=3} a user-defined prior is used, defined by \code{prior.dag}. It is given by an adjacency matrix (squared and same size as number of nodes) where entries ranging from zero to one give the user prior belief. An hyperparameter defining the global user belief in the prior is given by \code{prior.lambda}.
MCMC sampler comes with asymptotic statistical guarantees. Therefore it is highly advised to run multiple long enough chains. The burn-in phase length (i.e. throwing away first MCMC iterations) should be adequately chosen.
The argument \code{data.dists} must be a list with named arguments, one for each of the variables in \code{data.df}, where each entry is either \code{"poisson"}, \code{"binomial"}, or \code{"gaussian"}.
The parameter \code{heating} could improve convergence. It should be a real positive number. One is neutral. The larger, the more probable to accept any move.
}
\value{A list with an entry for the list of sampled DAGs, the list of scores, the acceptance probability, the method used for each MCMC jump, the rejection status for each MCMC jump, the total number of iterations the thinning, the length of burn-in phase, the named list of distribution per node, the heating parameter for each chain and a data.frame with the score of all chains. The returned object is of class mcmcabn.}
\author{Gilles Kratzer}
\references{
For the implementation of the function:
Kratzer, G., Furrer, R. "Is a single unique Bayesian network enough to accurately represent your data?". arXiv preprint arXiv:1902.06641.
For the new edge reversal:
Grzegorczyk, M., Husmeier, D. (2008). "Improving the structure MCMC sampler for Bayesian networks by introducing a new edge reversal move", Machine Learning, vol. 71(2-3), 265.
For the Markov Blanket resampling move:
Su, C., Borsuk, M. E. (2016). "Improving structure MCMC for Bayesian networks through Markov blanket resampling", The Journal of Machine Learning Research, vol. 17(1), 4042-4061.
For the Coupled Heated MCMC algorithm:
Geyer, C. J. (1991). Markov chain Monte Carlo maximum likelihood.
For the Koivisto prior:
Koivisto, M. V. (2004). Exact Structure Discovery in Bayesian Networks, Journal of Machine Learning Research, vol 5, 549-573.
For the user defined prior:
Werhli, A. V., Husmeier, D. (2007). "Reconstructing gene regulatory networks with Bayesian networks by combining expression data with multiple sources of prior knowledge". Statistical Applications in Genetics and Molecular Biology, 6 (Article 15).
Imoto, S., Higuchi, T., Goto, T., Tashiro, K., Kuhara, S., Miyano, S. (2003). Using Bayesian networks for estimating gene networks from microarrays and biological knowledge. In Proceedings of the European Conference on Computational Biology.
For the asia dataset:
Scutari, M. (2010). Learning Bayesian Networks with the bnlearn R Package. Journal of Statistical Software, 35(3), 1-22. doi:http://dx.doi.org/10.18637/jss.v035.i03.
}
\examples{
## Example from the asia dataset from Lauritzen and Spiegelhalter (1988)
## provided by Scutari (2010)
# The number of MCMC run is deliberately chosen too small (computing time)
# no thinning (usually not recommended)
# no burn-in (usually not recommended,
# even if not supported by any theoretical arguments)
data("mcmc_run_asia")
# run: 0.03 REV, 0.03 MBR, 0.94 MC3 MCMC jumps and 4 chains
# with a random DAG as starting point
mcmc.out.asia.small1 <- CoupledHeatedmcmcabn(score.cache = bsc.compute.asia,
score = "mlik",
data.dists = dist.asia,
max.parents = 2,
mcmc.scheme = c(150,0,0),
seed = 5416,
verbose = FALSE,
start.dag = "random",
prob.rev = 0.03,
prob.mbr = 0.03,
prior.choice = 2,heating = 0.,n.chains = 4)
summary(mcmc.out.asia.small1)
# compared to the mcmcabn() function
mcmc.out.asia.small2 <- mcmcabn(score.cache = bsc.compute.asia,
score = "mlik",
data.dists = dist.asia,
max.parents = 2,
mcmc.scheme = c(150,0,0),
seed = 5416,
verbose = FALSE,
start.dag = "random",
prob.rev = 0.03,
prob.mbr = 0.03,
prior.choice = 2,heating = 0.25)
summary(mcmc.out.asia.small2)
}
......@@ -489,6 +489,76 @@ expect_equal(query(mcmcabn = mcmc.out.gaussian.test)[2,3],query(mcmcabn = mcmc.o
expect_equal(query(mcmcabn = mcmc.out.gaussian.test)[6,1],query(mcmcabn = mcmc.out.gaussian.test,formula = ~G|A))
})
test_that("mcmcabnHeated",{
data(asia)
dist.asia <- list(Asia = "binomial",
Smoking = "binomial",
Tuberculosis = "binomial",
LungCancer = "binomial",
Bronchitis = "binomial",
Either = "binomial",
XRay = "binomial",
Dyspnea = "binomial")
colnames(asia) <- c("Asia","Smoking", "Tuberculosis", "LungCancer", "Bronchitis", "Either", "XRay", "Dyspnea")
bsc.compute.asia <- buildscorecache(data.df = asia,
data.dists = dist.asia,
max.parents = 2)
mcmc.out.asia <- CoupledHeatedmcmcabn(score.cache = bsc.compute.asia,
score = "mlik",
data.dists = dist.asia,
max.parents = 2,
mcmc.scheme = c(2000,0,0),
seed = 9,
verbose = FALSE,
start.dag = "random",
prob.rev = 0.1,
prob.mbr = 0.1,
prior.choice = 1,heating = 5,n.chains = 4)
#maximum scoring network using exact search (not MCMC based)
dag <- mostprobable(score.cache = bsc.compute.asia)
expect_equal(max(mcmc.out.asia$score.coupled),fitabn(object = dag,data.df <