Nothing
MBR <- function(n.var, dag.tmp, retain, ban, max.parents, sc, score.cache, score, verbose, heating) {
## scores
score.G <- 0
score.G.prime <- 0
rejection <- 1
A <- NULL
# store number of edges
n.edges <- sum(dag.tmp)
# Randomly select a node Xi in current graph G, withprobability N^−1
selected.node <- sample(x = 1:n.var, size = 1)
# store current parent set
parent.set <- dag.tmp[selected.node, ,drop=FALSE]
# remove parent set + co parent
dag.G.0 <- dag.tmp
# 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, ,drop=FALSE]
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]
# sample new parent set
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]
}
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
## new co parent set
if (sum(dag.G.1[, selected.node]) != 0) {
score.G <- numeric(sum(dag.G.1[, selected.node]))
# randomization of children
order.child <- sample(1:sum(dag.G.1[, selected.node]))
for (j in order.child) {
# descendant of node j
descendents.j <- descendents(nodes = j, dag = dag.G.1)
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == j, , drop=FALSE]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1,,drop=FALSE]
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendents.j,drop=FALSE] == 0) == length(descendents.j),, drop=FALSE]
# sample a new parent set
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]
}
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.MBR <- dag.G.1
############################ complementary inverse move
if(FALSE){
# parent set
parent.set.G.1 <- dag.G.1[selected.node, ]
# delete co-parent
tmp <- which(x = dag.G.1[, selected.node] == 1, arr.ind = TRUE)
if (!is.integer0(tmp)) {
dag.G.1[tmp[1], -selected.node] <- 0
}
# delete parents
dag.G.1[selected.node, ] <- 0
# from score take out parent + descendant
sc.tmp <- sc[score.cache$children == selected.node, ,drop=FALSE]
if (sum(parent.set.G.1) != 0) {
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)) {
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, ,drop=FALSE]
}
} else {
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,drop=FALSE] == 0) == length(descendents.i), ,drop=FALSE]
}
}
}
# store score
if (is.null(ncol(sc.tmp))) {
z.star.i.G.prime.0 <- sc.tmp[length(sc.tmp)]
} else {
z.star.i.G.prime.0 <- logSumExp(sc.tmp[, ncol(sc.tmp)])
}
# equivalent to G1
dag.G.1[selected.node, ] <- parent.set
if (sum(dag.G.1[, selected.node]) != 0) {
score.G.prime <- numeric(sum(dag.G.1[, selected.node]))
for (j in 1:sum(dag.G.1[, selected.node])) {
# 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, ,drop=FALSE]
sc.tmp <- sc.tmp[sc.tmp[, selected.node] == 1, ,drop=FALSE]
# remove descendents of j
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,drop=FALSE] == 0) == length(descendents.j), ]
}
}
}
dag.G.1[j, ] <- dag.tmp[j, ]
if (is.null(ncol(sc.tmp))) {
score.G.prime[j] <- logSumExp(sc.tmp[length(sc.tmp)])
} else {
score.G.prime[j] <- logSumExp(sc.tmp[, ncol(sc.tmp)])
}
}
}
}
}
}
############################## Acceptance probability
n.edges.tilde <- sum(dag.MBR)
#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))
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)
if(is.nan(A)){A <- 0}
if (runif(1)<exponent(A, heating)) {
rejection <- 0
dag.tmp <- dag.MBR
score <- score.dag(dag = dag.MBR,bsc.score = score.cache,sc = sc)
}
if (is.null(A)) {
A <- 0
}
############################## Return
return(list(dag.tmp = dag.tmp, score = score, alpha = A, rejection = rejection))
} #EOF
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.