Commit 33dbbeb5 authored by Gilles Kratzer's avatar Gilles Kratzer
Browse files

update

parent 469aa08d
Pipeline #2425 passed with stage
in 3 seconds
......@@ -14,6 +14,7 @@ print.mcmcabn
importFrom(cowplot, ggdraw, axis_canvas, insert_yaxis_grob)
importFrom(graphics, plot)
importFrom(stats, as.formula, dist, quantile, sd)
importFrom(matrixStats, logSumExp)
#methods for class mcmcabn
S3method(plot, mcmcabn)
......
MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
MBR <- function(n.var, dag.tmp, retain, ban, max.parents, sc, score.cache, score, verbose) {
## scores
score.G <- 0
......@@ -16,68 +16,40 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
# store current parent set
parent.set <- dag.tmp[selected.node, ]
parent.set <- dag.tmp[selected.node, ,drop=FALSE]
# remove parent set + co parent
dag.G.0 <- dag.tmp
tmp <- which(x = dag.G.0[, selected.node] == 1, arr.ind = TRUE)
if (!is.integer0(tmp)) {
dag.G.0[tmp[1], -selected.node] <- 0
}
dag.G.0[selected.node, ] <- 0
# descendant of node i
descendents.i <- descendents(nodes = selected.node, dag = dag.G.0)
#delete parent of children of node i
dag.G.0[descendents.i, -selected.node] <- 0
dag.G.0[selected.node, ] <- 0
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == selected.node, ]
if (sum(parent.set) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, c(as.logical(parent.set), FALSE)] == 0, ]
}
if (sum(parent.set) > 1) {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, c(as.logical(parent.set), FALSE)]) == sum(parent.set), ]
}
sc.tmp <- sc[score.cache$children == selected.node, ,drop=FALSE]
# remove descendents of i
if (!is.null(descendents.i)) {
if (length(descendents.i) == 1) {
if (is.null(dim(sc.tmp))) {
sc.tmp <- sc.tmp[sc.tmp[descendents.i] == 0]
} else {
sc.tmp <- sc.tmp[sc.tmp[, descendents.i] == 0, ]
}
sc.tmp <- sc.tmp[rowSums(sc.tmp[, c(as.logical(parent.set), FALSE),drop=FALSE]) == sum(parent.set), ,drop=FALSE]
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.i, drop=FALSE] == 0) == length(descendents.i), ,drop=FALSE]
} else {
# sample new parent set
if (is.null(dim(sc.tmp))) {
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.i] == 0) == length(descendents.i)]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.i] == 0) == length(descendents.i), ]
}
if(nrow(sc.tmp)==0){new.parent.i <- matrix(data = 0,nrow = 1,ncol = n.var+1)}
if(nrow(sc.tmp)==1){new.parent.i <- sc.tmp
}else{
new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] - sum(sc.tmp[, ncol(sc.tmp)]))), ,drop=FALSE]
}
}
# sample new parent set
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp))
{
if (is.vector(sc.tmp)) {
new.parent.i <- sc.tmp
new.parent.i <- new.parent.i[-length(new.parent.i)]
# store partition function
z.i.G.0 <- sc.tmp[length(sc.tmp)]
} else {
# sample a new parent set
new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
new.parent.i <- new.parent.i[-length(new.parent.i)]
z.i.G.0 <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
# dag G1
new.parent.i <- new.parent.i[,-length(new.parent.i)]
z.i.G.0 <- (logSumExp(sc.tmp[, ncol(sc.tmp)]))
dag.G.1 <- dag.G.0
dag.G.1[selected.node, ] <- new.parent.i
......@@ -92,57 +64,26 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
# descendant of node j
descendents.j <- descendents(nodes = j, dag = dag.G.1)
if (!is.character(descendents.j)) {
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == j, ]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1, ]
# print(class(sc.tmp)) remove descendents of i
if (is.null(ncol(sc.tmp))) {
if (!is.null(descendents.j)) {
if (length(descendents.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[descendents.j] == 0]
} else {
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.j] == 0) == length(descendents.j)]
}
}
} else {
if (!is.null(descendents.j)) {
if (length(descendents.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.j] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.j] == 0) == length(descendents.j), ]
}
}
}
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
if (!identical(unname(sc.tmp[length(sc.tmp)]), numeric(0))) {
sc.tmp <- sc[score.cache$children == j, , drop=FALSE]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1,,drop=FALSE]
if (is.vector(sc.tmp)) {
new.parent.j <- sc.tmp
new.parent.j <- new.parent.j[-length(new.parent.j)]
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.j,drop=FALSE] == 0) == length(descendents.j),, drop=FALSE]
# store partition function print(sc.tmp[length(sc.tmp)])
# sample a new parent set
score.G[j] <- sc.tmp[length(sc.tmp)]
if(nrow(sc.tmp)==0){new.parent.j <- matrix(data = 0,nrow = 1,ncol = n.var+1)}
if(nrow(sc.tmp)==1){new.parent.i <- sc.tmp}
if(nrow(sc.tmp)>1){
new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] - sum(sc.tmp[, ncol(sc.tmp)]))), ,drop=FALSE]
}
} else {
# sample a new parent set
new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
new.parent.j <- new.parent.j[-length(new.parent.j)]
new.parent.j <- new.parent.j[-length(new.parent.j)]
score.G[j] <- sum(sc.tmp[, ncol(sc.tmp)])
}
score.G[j] <- logSumExp(sc.tmp[, ncol(sc.tmp)])
dag.G.1[j, ] <- new.parent.j
}
}
}
}
}
......@@ -166,13 +107,10 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == selected.node, ]
sc.tmp <- sc[score.cache$children == selected.node, ,drop=FALSE]
if (sum(parent.set.G.1) != 0) {
if (sum(parent.set.G.1) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, c(as.logical(parent.set.G.1), FALSE)] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, c(as.logical(parent.set.G.1), FALSE)]) == sum(parent.set.G.1), ]
}
sc.tmp <- sc.tmp[rowSums(sc.tmp[, c(as.logical(parent.set.G.1), FALSE),drop=FALSE]) == sum(parent.set.G.1), ,drop=FALSE]
# remove descendents of i
if (!is.null(descendents.i)) {
......@@ -180,7 +118,7 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
if (is.null(dim(sc.tmp))) {
sc.tmp <- sc.tmp[sc.tmp[descendents.i] == 0]
} else {
sc.tmp <- sc.tmp[sc.tmp[, descendents.i] == 0, ]
sc.tmp <- sc.tmp[sc.tmp[, descendents.i] == 0, ,drop=FALSE]
}
} else {
......@@ -188,7 +126,7 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
if (is.null(dim(sc.tmp))) {
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.i] == 0) == length(descendents.i)]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.i] == 0) == length(descendents.i), ]
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.i,drop=FALSE] == 0) == length(descendents.i), ,drop=FALSE]
}
}
}
......@@ -198,7 +136,7 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
z.star.i.G.prime.0 <- sc.tmp[length(sc.tmp)]
} else {
z.star.i.G.prime.0 <- sum(sc.tmp[, ncol(sc.tmp)])
z.star.i.G.prime.0 <- logSumExp(sc.tmp[, ncol(sc.tmp)])
}
# equivalent to G1
......@@ -214,8 +152,8 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
descendents.j <- descendents(nodes = j, dag = dag.G.1)
if (!is.character(descendents.j)) {
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == j, ]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1, ]
sc.tmp <- sc[score.cache$children == j, ,drop=FALSE]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1, ,drop=FALSE]
# remove descendents of j
if (is.null(ncol(sc.tmp))) {
......@@ -232,7 +170,7 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
if (length(descendents.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.j] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.j] == 0) == length(descendents.j), ]
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.j,drop=FALSE] == 0) == length(descendents.j), ]
}
}
}
......@@ -240,9 +178,9 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
dag.G.1[j, ] <- dag.tmp[j, ]
if (is.null(ncol(sc.tmp))) {
score.G.prime[j] <- sum(sc.tmp[length(sc.tmp)])
score.G.prime[j] <- logSumExp(sc.tmp[length(sc.tmp)])
} else {
score.G.prime[j] <- sum(sc.tmp[, ncol(sc.tmp)])
score.G.prime[j] <- logSumExp(sc.tmp[, ncol(sc.tmp)])
}
}
......@@ -251,30 +189,22 @@ MBR <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
############################## Acceptance probability
score.MBR <- score.dag(dag = dag.MBR,bsc.score = score.cache,sc = sc)
#score.A <- exp(z.i.G.0 - z.star.i.G.prime.0 + sum(score.G) - sum(score.G.prime))
score.A <- exp(- score + score.MBR)
n.edges.tilde <- sum(dag.G.1)
score.A <- ( z.i.G.0 / z.star.i.G.prime.0 * sum(is.finite(score.G)) / sum(is.finite(score.G.prime)))
#score.A <- n.edges/n.edges.tilde * ((z.star.x.i.M.dot) / (z.star.x.j.M.dot) * (z.x.j.M.cross) / (z.x.i.M.tilde.cross))
if(is.nan(score.A)){score.A <- 0}
if((score.A)<0){score.A <- 0}
A <- min(1, score.A)
if (rbinom(n = 1, size = 1, prob = A) == 1) {
rejection <- 0
dag.tmp <- dag.MBR
#score.MBR <- 0
# for (a in 1:n.var) {
# sc.tmp <- sc[score.cache$children == a, ]
# score.MBR <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
# ])))), n.var + 1]), score.MBR)
# }
score <- score.MBR
score <- score.dag(dag = dag.MBR,bsc.score = score.cache,sc = sc)
}
}
} #EOIF
if (is.null(A)) {
A <- 0
}
......
REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
REV <- function(n.var, dag.tmp, retain, ban, max.parents, sc, score.cache, score, verbose) {
rejection <- 1
A <- 0
......@@ -9,18 +9,18 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
# randomly select one
if (sum(dag.tmp) != 0)
{
if (sum(dag.tmp) != 0){
selected.edge <- which(x = dag.tmp == 1, arr.ind = TRUE)[sample(x = 1:n.edges, size = 1), ]
#i->j
selected.edge <- which(x = ((dag.tmp*(1-retain)*(1-t(ban))) == 1), arr.ind = TRUE)[sample(x = 1:n.edges, size = 1), ,drop=FALSE]
# store current parent set
# store current parent set (j)
parent.set <- dag.tmp[selected.edge[1], ,drop=FALSE]
parent.set <- dag.tmp[selected.edge[1], ]
# remove parent set
# remove parent set (i and j)
dag.M.dot <- dag.tmp
dag.M.dot[selected.edge[1], ] <- 0
dag.M.dot[selected.edge[1],] <- 0
dag.M.dot[selected.edge[2],] <- 0
# store descendents j row i col
......@@ -29,47 +29,20 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
# mark all parent sets of i score that do not include j
sc.tmp <- sc[score.cache$children == selected.edge[2], ]
sc.tmp <- sc.tmp[sc.tmp[, selected.edge[1]] == 1, ]
# remove descendents of i
if (is.null(ncol(sc.tmp))) {
if (!is.null(descendents.M.dot.i)) {
if (length(descendents.M.dot.i) == 1) {
sc.tmp <- sc.tmp[sc.tmp[descendents.M.dot.i] == 0]
} else {
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.M.dot.i] == 0) == length(descendents.M.dot.i)]
}
}
} else {
if (!is.null(descendents.M.dot.i)) {
if (length(descendents.M.dot.i) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.M.dot.i] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.i] == 0) == length(descendents.M.dot.i), ]
}
}
}
sc.tmp <- sc[score.cache$children == selected.edge[2], ,drop=FALSE]
sc.tmp <- sc.tmp[sc.tmp[, selected.edge[1]] == 1, ,drop=FALSE]
# sample new parent set
#if (!is.null(descendents.M.dot.i)) {
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.i,drop=FALSE] == 0) == length(descendents.M.dot.i), ,drop=FALSE]
#}
if (is.vector(sc.tmp)) {
new.parent.i <- sc.tmp
new.parent.i <- new.parent.i[-length(new.parent.i)]
# store partition function
z.star.x.i.M.dot <- sc.tmp[length(sc.tmp)]
} else {
# !!!exp missing
new.parent.i <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
sum(sc.tmp[, ncol(sc.tmp)]))), ,drop=FALSE]
new.parent.i <- new.parent.i[-length(new.parent.i)]
# store partition function
z.star.x.i.M.dot <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
z.star.x.i.M.dot <- (logSumExp(sc.tmp[, ncol(sc.tmp)]))
# M direct sum
dag.M.cross <- dag.M.dot
......@@ -78,107 +51,59 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
# descendant of i
descendents.M.cross.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.cross)
descendents.M.cross.j <- descendents(nodes = unname(selected.edge[1]), dag = dag.M.cross)
#descendents.M.j <- descendents(nodes = unname(selected.edge[1]), dag = dag.tmp)
# score node j not descendant of M
sc.tmp <- sc[score.cache$children == selected.edge[1], ]
sc.tmp <- sc[score.cache$children == selected.edge[1], ,drop=FALSE]
# remove descendents of i
if (!is.null(descendents.M.cross.j)) {
if (length(descendents.M.cross.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.M.cross.j] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.cross.j] == 0) == length(descendents.M.cross.j),
]
}
}
#if (!is.null(descendents.M.cross.j)) {
#sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.j,drop=FALSE] == 0) == length(descendents.M.j),, drop=FALSE]
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.cross.j,drop=FALSE] == 0) == length(descendents.M.cross.j),, drop=FALSE]
#}
# sample new parent set
if (nrow(sc.tmp) > 0 || is.vector(sc.tmp)) {
if (is.vector(sc.tmp)) {
new.parent.j <- sc.tmp
# print(new.parent.i)
new.parent.j <- new.parent.j[-length(new.parent.j)]
# print(new.parent.i) store partition function
z.x.j.M.cross <- sc.tmp[length(sc.tmp)]
} else {
## !!! exp missing
new.parent.j <- sc.tmp[sample(x = 1:nrow(sc.tmp), size = 1, prob = range01(sc.tmp[, ncol(sc.tmp)] -
sum(sc.tmp[, ncol(sc.tmp)]))), ]
sum(sc.tmp[, ncol(sc.tmp)]))), ,drop=FALSE]
new.parent.j <- new.parent.j[-length(new.parent.j)]
# store partition function
z.x.j.M.cross <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
z.x.j.M.cross <- (logSumExp(sc.tmp[, ncol(sc.tmp)]))
# M tilde
dag.M.tilde <- dag.M.cross
dag.M.tilde[selected.edge[1], ] <- new.parent.j
dag.M.tilde[selected.edge[1],] <- new.parent.j
n.edges.tilde <- sum(dag.M.tilde)
############################## computing acceptance proba
############################## computing acceptance probability ##############################################
# score node j that do not include i
sc.tmp <- sc[score.cache$children == selected.edge[1], ]
sc.tmp <- sc.tmp[sc.tmp[, selected.edge[2]] == 1, ]
# remove descendents of i
if (is.vector(sc.tmp)) {
if (!is.null(descendents.M.dot.j)) {
if (length(descendents.M.dot.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[descendents.M.dot.j] == 0] #
} else {
# print(class(sc.tmp)) print(descendents.M.dot.j)
sc.tmp <- sc.tmp[sum(sc.tmp[descendents.M.dot.j] == 0) == length(descendents.M.dot.j)]
}
}
} else {
if (!is.null(descendents.M.dot.j)) {
if (length(descendents.M.dot.j) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.M.dot.j] == 0, ]
} else {
# print(class(sc.tmp)) print(descendents.M.dot.j)
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.j] == 0) == length(descendents.M.dot.j),
]
}
}
}
if (is.vector(sc.tmp)) {
z.star.x.j.M.dot <- sc.tmp[length(sc.tmp)]
} else {
z.star.x.j.M.dot <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
sc.tmp <- sc[score.cache$children == selected.edge[1], ,drop=FALSE]
sc.tmp <- sc.tmp[sc.tmp[, selected.edge[2]] == 1, ,drop=FALSE]
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.dot.j,drop=FALSE] == 0) == length(descendents.M.dot.j), ,drop=FALSE]
z.star.x.j.M.dot <- (logSumExp(sc.tmp[, ncol(sc.tmp)]))
dag.M.tilde.cross <- dag.M.dot
dag.M.tilde.cross[selected.edge[1], ] <- parent.set
dag.M.tilde.cross[selected.edge[1],] <- parent.set
# descendant
descendents.M.tilde.cross.i <- descendents(nodes = unname(selected.edge[2]), dag = dag.M.tilde.cross)
sc.tmp <- sc[score.cache$children == selected.edge[2], ]
sc.tmp <- sc[score.cache$children == selected.edge[2], ,drop=FALSE]
if (!is.null(descendents.M.tilde.cross.i)) {
if (length(descendents.M.tilde.cross.i) == 1) {
sc.tmp <- sc.tmp[sc.tmp[, descendents.M.tilde.cross.i] == 0, ]
} else {
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.tilde.cross.i] == 0) == length(descendents.M.tilde.cross.i),
]
sc.tmp <- sc.tmp[rowSums(sc.tmp[, descendents.M.tilde.cross.i,drop=FALSE] == 0) == length(descendents.M.tilde.cross.i),
,drop=FALSE]
}
}
z.x.i.M.tilde.cross <- (logSumExp(sc.tmp[, ncol(sc.tmp)]))
if (is.vector(sc.tmp)) {
z.x.i.M.tilde.cross <- sc.tmp[length(sc.tmp)]
} else {
z.x.i.M.tilde.cross <- (sum(sc.tmp[, ncol(sc.tmp)]))
}
############################## Acceptance probability
score.A <- n.edges/n.edges.tilde * exp(z.star.x.i.M.dot - z.star.x.j.M.dot + z.x.j.M.cross - z.x.i.M.tilde.cross)
score.A <- n.edges/n.edges.tilde * ((z.star.x.i.M.dot) / (z.star.x.j.M.dot) * (z.x.j.M.cross) / (z.x.i.M.tilde.cross))
if(is.nan(score.A)){score.A <- 0}
if((score.A)<0){score.A <- 0}
A <- min(1, score.A)
if (rbinom(n = 1, size = 1, prob = A) == 1) {
......@@ -187,15 +112,13 @@ REV <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, verbose) {
score.M.tilde <- 0
for (a in 1:n.var) {
sc.tmp <- sc[score.cache$children == a, ]
sc.tmp <- sc[score.cache$children == a, , drop=FALSE]
score.M.tilde <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
])))), n.var + 1]), score.M.tilde)
}
score <- score.M.tilde
}
}
}
} #EOFif
}#eoif
############################## Return
......
mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choice, prior.lambda, prior.dag, verbose) {
mc3 <- function(n.var, dag.tmp, retain, ban, max.parents, sc, score.cache, score, prior.choice, prior.lambda, prior.dag, verbose) {
## construction of neighbours list
neighbours.list <- NULL
......@@ -20,7 +20,7 @@ mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choic
## Removal of an arc
for (a in 1:n.var) {
for (b in 1:n.var) {
if (dag.tmp[a, b] == 1) {
if (dag.tmp[a, b] == 1 && retain[a, b]!=1) {
tmp <- dag.tmp
tmp[a, b] <- 0
neighbours.list[[length(neighbours.list) + 1]] <- tmp
......@@ -32,7 +32,7 @@ mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choic
## Addition of an arc
for (a in 1:n.var) {
for (b in 1:n.var) {
if ((dag.tmp[a, b] == 0 & dag.tmp[b, a] == 0) & sum(dag.tmp[a, ]) <= (max.parents - 1)) {
if(((dag.tmp[a, b] == 0 & dag.tmp[b, a] == 0) & sum(dag.tmp[a, ]) <= (max.parents - 1)) && ban[a,b]!=1) {
tmp <- dag.tmp
# print(tmp)
tmp[a, b] <- 1
......@@ -49,7 +49,7 @@ mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choic
## Reversal of an arc
for (a in 1:n.var) {
for (b in 1:n.var) {
if (dag.tmp[a, b] == 1 & sum(dag.tmp[b, ]) <= (max.parents - 1)) {
if ((dag.tmp[a, b] == 1 & sum(dag.tmp[b, ]) <= (max.parents - 1)) && retain[a,b] && ban[b,a]!=1) {
tmp <- dag.tmp
tmp[a, b] <- 0
tmp[b, a] <- 1
......@@ -133,7 +133,7 @@ mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choic
prior.Gprime <- prod(1/choose(n = (n.var - 1), k = sum(dag.gprime[a, ])), prior.Gprime)
}
sc.tmp <- sc[score.cache$children == a, ]
sc.tmp <- sc[score.cache$children == a,,drop = FALSE]
score.G <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.tmp[a,
])))), n.var + 1]), score.G)
score.Gprime <- sum(min(sc.tmp[(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.gprime[a,
......@@ -142,7 +142,7 @@ mc3 <- function(n.var, dag.tmp, max.parents, sc, score.cache, score, prior.choic
# 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)
if(!is.numeric(alpha)) alpha <- 0
if(!is.numeric(alpha) | is.nan(alpha)) alpha <- 0
......
......@@ -42,6 +42,16 @@ mcmcabn <- function(score.cache = NULL, score = "mlik", data.dists = NULL, max.p
## 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)
...