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

updates

parent 304162e0
MBR<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, verbose){
##scores
score.G <- 0
score.G.prime <- 0
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,]
#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
descendants.i <- descendants(nodes = selected.node,dag = dag.G.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),]
}
# remove descendants of i
if(!is.null(descendants.i)){
if(length(descendants.i)==1){
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sc.tmp[descendants.i]==0]
}else{
sc.tmp <- sc.tmp[sc.tmp[,descendants.i]==0,]
}
}else{
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sum(sc.tmp[descendants.i]==0)==length(descendants.i)]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.i]==0)==length(descendants.i),]}
}
}
#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
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
descendants.j <- descendants(nodes = j,dag = dag.G.1)
if(!is.character(descendants.j)){
#from score take out parent + descendant
sc.tmp <- sc[score.cache$children==j,]
sc.tmp <- sc.tmp[sc.tmp[,selected.node]==1,]
# remove descendants of i
if(!is.null(descendants.j)){
if(length(descendants.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.j]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.j]==0)==length(descendants.j),]}
}
if(nrow(sc.tmp)>0 || is.vector(sc.tmp)){
if(is.vector(sc.tmp)){
new.parent.j <- sc.tmp
new.parent.j<-new.parent.j[-length(new.parent.j)]
#store partition function
score.G[j] <- sc.tmp[length(sc.tmp)]
}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)]
score.G[j] <- sum(sc.tmp[,ncol(sc.tmp)])
}
dag.G.1[j,] <- new.parent.j
}
}}}
dag.MBR <- dag.G.1
############################
##complementary inverse move
############################
#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,]
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),]
}
# remove descendants of i
if(!is.null(descendants.i)){
if(length(descendants.i)==1){
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sc.tmp[descendants.i]==0]
}else{
sc.tmp <- sc.tmp[sc.tmp[,descendants.i]==0,]
}
}else{
if(is.null(dim(sc.tmp))){
sc.tmp <- sc.tmp[sum(sc.tmp[descendants.i]==0)==length(descendants.i)]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.i]==0)==length(descendants.i),]}
}
}
#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 <- sum(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
descendants.j <- descendants(nodes = j,dag = dag.G.1)
if(!is.character(descendants.j)){
#from score take out parent + descendant
sc.tmp <- sc[score.cache$children==j,]
sc.tmp <- sc.tmp[sc.tmp[,selected.node]==1,]
# remove descendants of j
if(!is.null(descendants.j)){
if(length(descendants.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.j]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.j]==0)==length(descendants.j),]}
}
dag.G.1[j,] <- dag.tmp[j,]
score.G.prime[j] <- sum(sc.tmp[,ncol(sc.tmp)])
}}}
##############################
## Acceptance probability
##############################
score.A <- exp(z.i.G.0-z.star.i.G.prime.0 +sum(score.G)-sum(score.G.prime))
A <- min(1,score.A)
if(rbinom(n = 1,size = 1,prob = A)==1){
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
}
}}#EOIF
if(is.null(A)){A<-0}
##############################
## Return
##############################
return(list(dag.tmp = dag.tmp, score = score, alpha=A))
}#EOF
......@@ -20,11 +20,27 @@ mcmc <- function(score.cache = NULL,
start.dag = NULL,
dag.retained = NULL,
dag.banned = NULL,
prob.rev = 0.15,
prob.rev = 0.05,
prob.mbr = 0.05,
prior = 2){
####################################################
##Tests
####################################################
if(is.null(score.cache))stop("A cache of score should be provided. YOu can produce it using the R package abn.")
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(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).")
##end of tests
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 <- rep(NA,length = mcmc.scheme[1]+1)
......@@ -90,27 +106,30 @@ mcmc <- function(score.cache = NULL,
while(j<=mcmc.scheme[3]){
####################################################
##REV move
####################################################
if(rbinom(n = 1,size = 1,prob = prob.rev)==1 & sum(dag.tmp)!=0){
if(verbose){print("REV move")}
out <- REV(n.var,(dag.tmp),max.parents,sc,score.cache,score,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
##method choice:
j<-j+1
}#EOF: REV move
else{
method.choice <- sample(x = c("mc3","REV","MBR"),size = 1,prob = c(prob.mc3,prob.rev,prob.mbr))
out <- mc3(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
switch (method.choice, "mc3" = {
out <- mc3(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
},
"REV" = {
if(verbose){print("REV move")}
out <- REV(n.var,(dag.tmp),max.parents,sc,score.cache,score,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
},
"MBR" = {
if(verbose){print("MBR move")}
out <- MBR(n.var,(dag.tmp),max.parents,sc,score.cache,score,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
})
j<-j+1
}
}
}#EOF
out.scores[1] <- score
......@@ -124,24 +143,30 @@ mcmc <- function(score.cache = NULL,
j <- 0
while(j<=mcmc.scheme[2]){
if(rbinom(n = 1,size = 1,prob = prob.rev)==1){
if(verbose){print("REV move")}
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),max.parents,sc,score.cache,score,prior,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
alpha <- out$alpha
},
"REV" = {
if(verbose){print("REV move")}
out <- REV(n.var,(dag.tmp),max.parents,sc,score.cache,score,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
alpha <- out$alpha
j <- j+1 ## while update
}#EOF: REV move
else{
out <- mc3(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior,verbose)
},
"MBR" = {
if(verbose){print("MBR move")}
out <- MBR(n.var,(dag.tmp),max.parents,sc,score.cache,score,verbose)
dag.tmp <- out$dag.tmp
score <- out$score
alpha <- out$alpha
j <- j+1 ## while update
}
})
j<-j+1
}#EOF mcmc search
out.dags[,,mcmc.search+1] <- dag.tmp
......
......@@ -136,21 +136,28 @@ is.integer0 <- function(x)
is.integer(x) && length(x) == 0L
}
ancestors<-function(nodes, dag){
if(!is.integer0(which(x = dag[nodes,]==1,arr.ind = TRUE))){
tmp <- which(x = dag[nodes,]==1,arr.ind = TRUE)
return(unname(c(tmp,ancestors(nodes = tmp,dag = dag))))
}else{
return(NULL)
}
}
##-------------------------------------------------------------------------
##Ancestor function
##-------------------------------------------------------------------------
#
# ancestors<-function(nodes, dag){
# diag(dag)<-0
# if(!is.integer0(which(x = dag[nodes,]==1,arr.ind = TRUE))){
# tmp <- which(x = dag[nodes,]==1,arr.ind = TRUE)
# return(unname(c(tmp,ancestors(nodes = tmp,dag = dag))))
# }else{
# return(NULL)
# }
# }
##-------------------------------------------------------------------------
##Descendant function
##-------------------------------------------------------------------------
descendants<-function(nodes, dag){
diag(dag)<-0
if(!identical(topoSortMAT(amat = dag,index = FALSE),character(0))){
if(!is.integer0(which(x = dag[,nodes]==1,arr.ind = TRUE)) && diag(dag)==0){
if(length(nodes)==1){
tmp <- unname(which(x = dag[,nodes]==1,arr.ind = TRUE))}else{
......@@ -160,6 +167,10 @@ descendants<-function(nodes, dag){
}else{
return(NULL)
}
}
}
return("Not a DAG")
}
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
'%!in%' <- function(x,y)!('%in%'(x,y))
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment