Commit a90e9cc8 authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

tests + priors

parent be36125e
Pipeline #1578 passed with stage
in 2 seconds
......@@ -20,4 +20,4 @@ Suggests: bnlearn, knitr, rmarkdown, ggdag
RoxygenNote: 6.1.1
VignetteBuilder: knitr
URL: https://www.math.uzh.ch/pages/mcmcabn/
BugReports: https://git.math.uzh.ch/gkratz/mcmcabn/issues
\ No newline at end of file
BugReports: https://git.math.uzh.ch/gkratz/mcmcabn/issues
REV<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, verbose){
rejection <- 1
A <- 0
#stor number of edges
......@@ -8,6 +9,8 @@ n.edges <- sum(dag.tmp)
#randomly select one
if(sum(dag.tmp)!=0){
selected.edge <- which(x = dag.tmp==1,arr.ind = TRUE)[sample(x = 1:n.edges,size = 1),]
#store current parent set
......@@ -183,8 +186,8 @@ if(rbinom(n = 1,size = 1,prob = A)==1){
}
score <- score.M.tilde
}
}}#EOFif
if(is.null(A)){A <- 0}
}}}#EOFif
##############################
## Return
......
......@@ -114,18 +114,18 @@ mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior.choice,
##prior computation
prior.G <- prior.Gprime <- 1
score.G <- score.Gprime <- 0
for(a in 1:n.var){
#user prior
if(prior.choice==3){
prior.G <- exp(-prior.lambda*sum(abs(dag.tmp-prior.dag)))
prior.Gprime <-exp(-prior.lambda*sum(abs(dag.gprime-prior.dag)))
}
for(a in 1:n.var){
#Koivisto priors
if(prior.choice==2){
prior.G <- sum(1/choose(n = (n.var-1),k = sum(dag.tmp[a,])),prior.G)
prior.Gprime <- sum(1/choose(n = (n.var-1),k = sum(dag.gprime[a,])),prior.Gprime)
}
#user prior
if(prior.choice==3){
prior.G <- sum(exp(-prior.lambda*abs(dag.tmp[a,]-prior.dag[a,])),prior.G)
prior.Gprime <-sum(exp(-prior.lambda*abs(dag.gprime[a,]-prior.dag[a,])),prior.Gprime)
prior.G <- prod(1/choose(n = (n.var-1),k = sum(dag.tmp[a,])),prior.G)
prior.Gprime <- prod(1/choose(n = (n.var-1),k = sum(dag.gprime[a,])),prior.Gprime)
}
sc.tmp <- sc[score.cache$children==a,]
......@@ -133,12 +133,14 @@ mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior.choice,
score.Gprime <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.gprime[a,])))),n.var+1]),score.Gprime)
}
alpha <- min(exp((score.Gprime + log(prior.Gprime) + log(n.Gprime)) - ((score.G + log(prior.G)) + log(n.G))),1)
#alpha <- min(exp((score.Gprime + log(prior.Gprime) + log(n.Gprime)) - ((score.G + log(prior.G)) + log(n.G))),1)
alpha <- min(exp(score.Gprime - score.G)*(n.G/n.Gprime)*(prior.Gprime/prior.G),1)
score <- score.G
if(!is.null(dag.gprime) & rbinom(n = 1,size = 1,prob = alpha)==1){
if(!is.null(dag.gprime) && rbinom(n = 1,size = 1,prob = alpha)==1){
dag.tmp <- dag.gprime
score <- score.Gprime
rejection <- 0
......
......@@ -33,15 +33,20 @@ mcmcabn <- function(score.cache = NULL,
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(max(rowSums(score.cache$node.defn))<max.parents)stop("Check max.parents. It should be the same as the one used in abn::buildscorecache() R function")
if(length(mcmc.scheme)!=3)stop("An MCMC scheme have to be provided. It should be such that c(returned,thinned,burned).")
if(length(mcmc.scheme)!=3)stop("An MCMC scheme have to be provided. It should be such that c(returned,thinned,burned) made of non negative integers.")
if(!is.numeric(mcmc.scheme[1]) | !is.numeric(mcmc.scheme[2]) | !is.numeric(mcmc.scheme[3]) | mcmc.scheme[1]<0 | mcmc.scheme[2]<0 | mcmc.scheme[3]<0){stop("An MCMC scheme have to be provided. It should be such that c(returned,thinned,burned) made of non negative integers.")}
if(max.parents<1 || max.parents>length(data.dists))stop("max.parents makes no sense.")
if(is.numeric(prob.rev) && prob.rev>1 && prob.rev<0)stop("prob.rev should be a probability.")
if(is.numeric(prob.mbr) && prob.mbr>1 && prob.mbr<0)stop("prob.mbr should be a probability.")
if(is.matrix(start.dag) && dim(start.dag)!=c(length(data.dists),length(data.dists)))stop("start.dag should be a squared matrix with dimension equal to the number of variables.")
if(is.matrix(start.dag) && is.null(topoSortMAT((start.dag)))){stop("start.dag should be a named DAG.")}
if(length(data.dists)!=max(score.cache$children))stop("data.dists should be a named list of all variables used to build the cache of precomputed score.")
if(length(data.dists)==0 || is.null(names(data.dists))){stop("data.dists should be a named list of all variables used to build the cache of precomputed score.")}
if(is.matrix(prior.dag)){prior.choice <- 3}
if(prior.choice %!in% 1:3){stop("prior.choice should be either 1,2 or 3.")}
if(is.matrix(prior.dag) && is.null(prior.lambda)){prior.lambda <- 1}
if(is.matrix(prior.dag) && dim(prior.dag)==c(length(data.dists),length(data.dists)))stop("prior.dag should be a squared matrix with dimension equal to the number of variables.")
if(is.matrix(prior.dag) && dim(prior.dag)!=c(length(data.dists),length(data.dists)))stop("prior.dag should be a squared matrix with dimension equal to the number of variables.")
##end of tests
n.var <- length(data.dists)
......@@ -68,15 +73,20 @@ mcmcabn <- function(score.cache = NULL,
##start random matrix
if(is.null(start.dag)){start.dag <- "random"}
if(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)]
dag.tmp <- dag.tmp[order(name.dag),order(name.dag)]
}
if(is.matrix(start.dag)){dag.tmp <- start.dag;colnames(dag.tmp) <- rownames(dag.tmp)<-1:n.var}
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")){
......
......@@ -21,7 +21,7 @@ ___
# mcmcabn: a structural MCMC sampler for DAGs learned from observed systemic datasets
## Installation
## Quick start
You can install `mcmabn` with:
......
......@@ -112,7 +112,7 @@
<p><strong>What are the R package functionalities?</strong></p>
<p>The workhorse of the R package is the <code><a href="../reference/mcmc.html">mcmcabn()</a></code> function. It generates MCMC samples from the posterior distribution of the data supported DAGs. This object can be summarized with a comprehensive verbal description or a plot. From those samples one can either compute the most probable structure and plot the cumulative maximum score achieved or query it using <code><a href="../reference/query.html">query()</a></code> R function using formula statement. For example what is the probability have an arc between two nodes, one node being in the Markov blanket of another or what is the probability to <strong>not</strong> have an arc between two nodes? Alternatively, one can query the frequency of an arc in one direction relative to the opposite direction.</p>
<p><strong>What is the structure of this document?</strong></p>
<p>We first illustrate the package with a simple example. Then some more theoretical details are provided. For a full description of the technical details we refer to the original publication [XXX].</p>
<p>We first illustrate the package with a simple example. Then some more theoretical details are provided. For a full description of the technical details we refer to the original publication.</p>
<div id="simple-mcmcabn-example" class="section level1">
<h1 class="hasAnchor">
<a href="#simple-mcmcabn-example" class="anchor"></a>Simple <em>mcmcabn</em> example</h1>
......@@ -120,7 +120,7 @@
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="https://www.rdocumentation.org/packages/utils/topics/install.packages">install.packages</a></span>(<span class="st">"mcmcabn"</span>)</code></pre></div>
<p>Once installed, the <code>mcmcabn</code> R package can be loaded using:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/library">library</a></span>(mcmcabn)</code></pre></div>
<p>Let us start with a ranking example from the <code>bnlearn</code> R package [XXX]. The first example is about a small synthetic dataset from Lauritzen and Spiegelhalter (1988) about lung diseases (tuberculosis, lung cancer or bronchitis) and visits to Asia.</p>
<p>Let us start with a ranking example from the <code>bnlearn</code> R package from Scutari (2010). The first example is about a small synthetic dataset from Lauritzen and Spiegelhalter (1988) about lung diseases (tuberculosis, lung cancer or bronchitis) and visits to Asia.</p>
<p>First thing first, one need to precompute a cache of scores. We use the R package <code>abn</code> to do it.</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/library">library</a></span>(bnlearn) <span class="co">#for the dataset</span>
<span class="kw"><a href="https://www.rdocumentation.org/packages/base/topics/library">library</a></span>(abn) <span class="co">#to precompute the scores </span>
......@@ -312,7 +312,7 @@ dag &lt;-<span class="st"> </span><span class="kw"><a href="https://www.rdocumen
<div id="priors" class="section level3">
<h3 class="hasAnchor">
<a href="#priors" class="anchor"></a>Priors</h3>
<p>Three priors are implemented in mcmcabn R package. The parameter <code>prior.choice</code> determines the prior used within each individual node for a given choice of parent combination. In Koivisto and Sood (2004) p.554 a form of prior, called Koivisto 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 <span class="math inline">\(|G_n|\)</span> is proportional to 1/[N-1 choose |G|] where there are N nodes in total.</p>
<p>Three priors are implemented in <code>mcmcabn</code> R package. The parameter <code>prior.choice</code> determines the prior used within each individual node for a given choice of parent combination. In Koivisto and Sood (2004) p.554 a form of prior, called Koivisto 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 <span class="math inline">\(|G_n|\)</span> is proportional to 1/[N-1 choose |G|] where there are N nodes in total.</p>
<p>Explicitly, the Koivisto prior is given by</p>
<p><span class="math display">\[p(G) = \frac{1}{z}\prod_{n=1}^N {N-1 \choose |G_n|}^{-1}\]</span></p>
<p>As we are in a MCMC setting and we will compute an Hasting ratio, the normalizing constant <span class="math inline">\(z\)</span> will cancelled out. Then we will never explicitly compute it.</p>
......@@ -342,7 +342,8 @@ Grzegorczyk, M.Husmeier, D. “Improving the structure MCMC sampler for Bayesian
Su, C.Borsuk, M. E. “Improving structure MCMC for Bayesian networks through Markov blanket resampling”, The Journal of Machine Learning Research, vol. 17(1), pp. 4042-4061, 2016.
Koivisto, M. V. (2004). Exact Structure Discovery in Bayesian Networks, Journal of Machine Learning Research, vol 5, 549-573.
Werhli, A. V., &amp; 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., &amp; Miyano, S. (2003). Using Bayesian networks for estimating gene networks from microarrays and biological knowledge. In Proceedings of the European Conference on Computational Biology.</div>
Imoto, S., Higuchi, T., Goto, T., Tashiro, K., Kuhara, S., &amp; Miyano, S. (2003). Using Bayesian networks for estimating gene networks from microarrays and biological knowledge. In Proceedings of the European Conference on Computational Biology.
Scutari, M. (2010). Learning Bayesian Networks with the bnlearn R Package. Journal of Statistical Software, 35(3), 1 - 22. <a href="doi:http://dx.doi.org/10.18637/jss.v035.i03" class="uri">doi:http://dx.doi.org/10.18637/jss.v035.i03</a>.</div>
</div>
</div>
......
......@@ -96,9 +96,9 @@
<div id="mcmcabn-a-structural-mcmc-sampler-for-dags-learned-from-observed-systemic-datasets" class="section level1">
<div class="page-header"><h1 class="hasAnchor">
<a href="#mcmcabn-a-structural-mcmc-sampler-for-dags-learned-from-observed-systemic-datasets" class="anchor"></a>mcmcabn: a structural MCMC sampler for DAGs learned from observed systemic datasets</h1></div>
<div id="installation" class="section level2">
<div id="quick-start" class="section level2">
<h2 class="hasAnchor">
<a href="#installation" class="anchor"></a>Installation</h2>
<a href="#quick-start" class="anchor"></a>Quick start</h2>
<p>You can install <code>mcmabn</code> with:</p>
<div class="sourceCode"><pre class="sourceCode r"><code class="sourceCode r"><span class="kw"><a href="https://www.rdocumentation.org/packages/utils/topics/install.packages">install.packages</a></span>(<span class="st">"mcmcabn"</span>)</code></pre></div>
<p>The three main problems addressed by this R package are:</p>
......
model{
###-----------------------
###Binomial nodes
###-----------------------
A ~ dbern(p.A); #Binary response
logit(p.A) <- 0 + 1*C; #Logistic regression
###-----------------------
###Binomial nodes
###-----------------------
B ~ dbern(p.B); #Binary response
logit(p.B) <- 1 + 1*A + 1*E; #Logistic regression
###-----------------------
###Binomial nodes
###-----------------------
C ~ dbern(p.C); #Binary response
logit(p.C) <- 0; #Logistic regression
###-----------------------
###Binomial nodes
###-----------------------
D ~ dbern(p.D); #Binary response
logit(p.D) <- 0 + 1*E; #Logistic regression
###-----------------------
###Binomial nodes
###-----------------------
E ~ dbern(p.E); #Binary response
logit(p.E) <- 1 + 1*C; #Logistic regression
}
\ No newline at end of file
###############################################################################
## test.R ---
## Author : Gilles Kratzer
## Document created: 14/11/2018
## Document created: 25/02/2019
###############################################################################
##Purpose: Test the mcmcabn sofware
##Purpose: Test the mcmcabn software
Sys.setenv("R_TESTS" = "")
library(testthat)
library(abn)
library(bnlearn)
##~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
##General tests
......@@ -17,23 +19,341 @@ context("General tests")
test_that("mc3",{
#add arrow
dist <- list(A="binomial",
B="binomial",
C="binomial",
D="binomial",
E="binomial")
data.param <- matrix(data = c(0,0,0.9,0,0,0.9,0,0.9,0,0.9,0,0,0,0,0,0,0,0,0,0.9,0,0,0.9,0,0.9),nrow = 5L,ncol = 5L,byrow = T)
diag(data.param)<-0.5
colnames(data.param) <- rownames(data.param) <- names(dist)
out.sim.0 <- invisible(simulateabn(data.dists = dist,n.chains = 1,n.adapt = 1000,n.thin = 1,n.iter = 1000,data.param = 0.4*data.param, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 3,
mcmc.scheme = c(1000,0,0),
seed = 42,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.choice = 1)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitabn(dag.m = dag,data.df = out.sim.0,data.dists = dist)$mlik)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 3,
mcmc.scheme = c(5000,0,0),
seed = 56,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.choice = 2)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 2)
expect_equal(round(max(mc3.out$scores),digits = 0),round(fitabn(dag.m = dag,data.df = out.sim.0,data.dists = dist)$mlik,digits = 0))
data.param.eq <- matrix(data = 0,nrow = 5,ncol = 5)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 1,
mcmc.scheme = c(1000,0,0),
seed = 42,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.dag = data.param.eq,
prior.lambda = 5,
prior.choice = 3)
expect_false(table(apply(mc3.out$dags,3, sum))[1]<500)
data.param.eq <- matrix(data = 1,nrow = 5,ncol = 5)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 1,
mcmc.scheme = c(1000,0,0),
seed = 42,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.dag = data.param.eq,
prior.lambda = 5,
prior.choice = 3)
expect_false(table(apply(mc3.out$dags,3, sum))[1]>500)
})
test_that("REV",{
dist <- list(A="binomial",
B="binomial",
C="binomial",
D="binomial",
E="binomial")
data.param <- matrix(data = c(0,0.9,0,0.9,0,0,0,0,0,0,0.9,0.9,0,0,0,0,0,0,0.9,0,0,0,0,0.9,0),nrow = 5L,ncol = 5L,byrow = T)
diag(data.param)<-0.5
colnames(data.param) <- rownames(data.param) <- names(dist)
out.sim.0 <- invisible(simulateabn(data.dists = dist,n.chains = 1,n.adapt = 1000,n.thin = 1,n.iter = 1000,data.param = 0.4*data.param, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 3,
mcmc.scheme = c(100,0,0),
seed = 23213,
verbose = FALSE,
start.dag = "random",
prob.rev = 0.5,
prob.mbr = 0,
prior.choice = 1)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitabn(dag.m = dag,data.df = out.sim.0,data.dists = dist)$mlik)
expect_silent(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 3,
mcmc.scheme = c(100,0,0),
seed = 32132,
verbose = FALSE,
start.dag = "random",
prob.rev = 1,
prob.mbr = 0,
prior.choice = 1))
})
test_that("MBR",{
dist <- list(A="binomial",
B="binomial",
C="binomial",
D="binomial",
E="binomial")
data.param <- matrix(data = c(0,0,0.9,0.9,0,0,0.9,0,0,0.9,0,0,0,0,0,0,0,0,0,0.9,0,0,0,0,0),nrow = 5L,ncol = 5L,byrow = T)
diag(data.param)<-0.5
colnames(data.param) <- rownames(data.param) <- names(dist)
out.sim.0 <- invisible(simulateabn(data.dists = dist,n.chains = 1,n.adapt = 1000,n.thin = 1,n.iter = 1000,data.param = 0.4*data.param, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 3)
mc3.out <- mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 3,
mcmc.scheme = c(100,0,0),
seed = 23213,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0.1,
prior.choice = 2)
dag <- mostprobable(score.cache = bsc.compute.0,prior.choice = 1)
expect_equal(max(mc3.out$scores),fitabn(dag.m = dag,data.df = out.sim.0,data.dists = dist)$mlik)
expect_silent(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 3,
mcmc.scheme = c(1000,0,0),
seed = 32132,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 1,
prior.choice = 1))
})
test_that("mcmcabn",{
dist<-list(A="binomial",
B="binomial",
C="binomial",
D="binomial",
E="binomial")
})
data.param.0 <- matrix(data = c(0,0,1,0,0,1,1,0,0,1,0,0,0,0,0,0,0,0,0,1,0,0,1,0,1),nrow = 5L,ncol = 5L,byrow = T)
colnames(data.param.0) <- rownames(data.param.0) <- names(dist)
out.sim.0 <- invisible(simulateabn(data.dists = dist,n.chains = 1,n.adapt = 20,n.thin = 1,n.iter = 100,data.param = data.param.0, simulate = TRUE,seed = 132))
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 2)
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "blabla",
data.dists = dist,
max.parents = 1,
mcmc.scheme = c(5,0,0),
seed = 2343,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.choice = 1))
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = unname(dist),
max.parents = 1,
mcmc.scheme = c(5,0,0),
seed = 2343,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.choice = 1))
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = unname(dist),
max.parents = 0,
mcmc.scheme = c(5,0,0),
seed = 2343,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.choice = 1))
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 0,
mcmc.scheme = c(-1,0,0),
seed = 2343,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.choice = 1))
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 1,
mcmc.scheme = c(1,0,0),
seed = 2343,
verbose = FALSE,
start.dag = "random",
prob.rev = 0,
prob.mbr = 0,
prior.choice = 15))
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 1,
mcmc.scheme = c(1,0,0),
seed = 2343,
verbose = FALSE,
start.dag = "random",
prob.rev = -0.1,
prob.mbr = 0,
prior.choice = 1))
expect_error(mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 1,
mcmc.scheme = c(1,0,0),
seed = 2343,
verbose = FALSE,
start.dag = "random",
prob.rev = -0.1,
prob.mbr = 0,
prior.choice = 1))
##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 <- mcmcabn(score.cache = bsc.compute.asia,
score = "mlik",
data.dists = dist.asia,
max.parents = 2,
mcmc.scheme = c(1000,0,0),
seed = 42,
verbose = FALSE,
start.dag = "random",
prob.rev = 0.03,
prob.mbr = 0.03,
prior.choice = 2)
max(mcmc.out.asia$scores)
#maximum scoring network using exact search (not MCMC based)
dag <- mostprobable(score.cache = bsc.compute.asia)
fitabn(dag.m = dag,data.df = asia,data.dists = dist.asia)$mlik
data.param.var.0 <- matrix(data = 0,nrow = 5L,ncol = 5L,byrow = T)
diag(data.param.var.0)<-0.1
dag<-plotabn(dag.m = data.param.0,data.dists = dist,plot = FALSE)
colnames(dag)<-rownames(dag)<-1:5
out.sim.0 <- simulateabn(data.dists = dist,n.chains = 1,n.adapt = 20,n.thin = 1,n.iter = 100,data.param = data.param.0, simulate = TRUE,seed = 132)
bsc.compute.0 <- buildscorecache(data.df = out.sim.0, data.dists = dist, max.parents = 2)
mc3.out.0<-mcmcabn(score.cache = bsc.compute.0,
score = "mlik",
data.dists = dist,
max.parents = 2,
mcmc.scheme = c(1000,1,1000),
seed = 2343,
verbose = TRUE,
start.dag = "random",
#algo = "mc3",
prob.rev = 0,
prob.mbr = 0,
prior = 2)
})