Commit 304162e0 authored by Reinhard Furrer's avatar Reinhard Furrer
Browse files

Merge branch 'master' of git.math.uzh.ch:gkratz/mcmcabn

parents f22005cb fb5ca7a3
REV<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, prior, verbose){
REV<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, verbose){
#stor number of edges
......@@ -6,7 +6,7 @@ n.edges <- sum(dag.tmp)
#randomly select one
selected.edge <- which(x = dag.tmp==1,arr.ind = TRUE)[sample(x = 1:n.edges,1),]
selected.edge <- which(x = dag.tmp==1,arr.ind = TRUE)[sample(x = 1:n.edges,size = 1),]
#store current parent set
......@@ -16,183 +16,149 @@ parent.set <- dag.tmp[selected.edge[1],]
dag.M.dot <- dag.tmp
dag.M.dot[selected.edge[1],] <- 0
ancestors()
#store descendant
# j row
# i col
descendants.M.dot.i <- descendants(nodes = unname(selected.edge[2]),dag = dag.M.dot)
descendants.M.dot.j <- descendants(nodes = unname(selected.edge[1]),dag = dag.M.dot)
#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 descendants of i
if(!is.null(descendants.M.dot.i)){
if(length(descendants.M.dot.i)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.dot.i]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.M.dot.i]==0)==length(descendants.M.dot.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.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)]))),]
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)]))
}
#M direct sum
dag.M.cross <- dag.M.dot
dag.M.cross[selected.edge[2],] <- new.parent.i
#descendant of i
descendants.M.cross.i <- descendants(nodes = unname(selected.edge[2]),dag = dag.M.cross)
descendants.M.cross.j <- descendants(nodes = unname(selected.edge[1]),dag = dag.M.cross)
# score node j not descendant of M
sc.tmp <- sc[score.cache$children==selected.edge[1],]
# remove descendants of i
if(!is.null(descendants.M.cross.j)){
if(length(descendants.M.cross.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.cross.j]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.M.cross.j]==0)==length(descendants.M.cross.j),]}
}
#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)]))),]
new.parent.j <- new.parent.j[-length(new.parent.j)]
#store partition function
z.x.j.M.cross <- (sum(sc.tmp[,ncol(sc.tmp)]))
}
#M tilde
dag.M.tilde <- dag.M.cross
dag.M.tilde[selected.edge[1],] <- new.parent.j
n.edges.tilde <- sum(dag.M.tilde)
############################## computing acceptance proba
# 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 descendants of i
if(!is.null(descendants.M.dot.j)){
if(length(descendants.M.dot.j)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.dot.j]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.M.dot.j]==0)==length(descendants.M.dot.j),]}
}
##old code
##sum of logs
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)]))
}
sumlogs<-function(x){
xprime<-sort(x,decreasing = TRUE)
n<-length(xprime)
out1<-0
if(n>1){
for(i in 2:n){
out1<-out1+exp(1)**(xprime[i]-xprime[1])
}}
out<-xprime[1]+log(1+out1)
dag.M.tilde.cross <- dag.M.dot
dag.M.tilde.cross[selected.edge[1],] <- parent.set
return(out)
}
#descendant
descendants.M.tilde.cross.i <- descendants(nodes = unname(selected.edge[2]),dag = dag.M.tilde.cross)
alpha <- 0
indices <- NULL
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)){
tmp <- dag.tmp
#print(tmp)
tmp[a,b] <- 0
tmp[b,a] <- 1
if(!identical(topoSortMAT(amat = tmp,index = FALSE),character(0))){
#neighbours.list[[length(neighbours.list)+1]] <- tmp
indices[[length(indices)+1]] <- c(a,b)
}
}
}
sc.tmp <- sc[score.cache$children==selected.edge[2],]
if(!is.null(descendants.M.tilde.cross.i)){
if(length(descendants.M.tilde.cross.i)==1){
sc.tmp <- sc.tmp[sc.tmp[,descendants.M.tilde.cross.i]==0,]
}else{
sc.tmp <- sc.tmp[rowSums(sc.tmp[,descendants.M.tilde.cross.i]==0)==length(descendants.M.tilde.cross.i),]}
}
if(length(indices)!=0){
##step1:
n.G <- length(indices)
mcmc.move <- sample(x = 1:n.G,size = 1,replace = FALSE)
dag.gprime <- dag.tmp
#reverse computing
##Q1:
sc.tmp <- sc[score.cache$children==indices[[mcmc.move]][1],]
sc.tmp <- sc.tmp[sc.tmp[,indices[[mcmc.move]][2]]==1,]
#Q1.tmp <- sc.tmp[which(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.gprime[indices[[mcmc.move]][1],])))),ncol(sc.tmp)]
#Q.G1 <- Q1.tmp/sum(sc.tmp[,ncol(sc.tmp)])
Q.G1 <- sumlogs(sc.tmp[,ncol(sc.tmp)])
##Q2:
sc.tmp <- sc[score.cache$children==indices[[mcmc.move]][2],]
#Q2.tmp <- sc.tmp[which(apply(sc.tmp, 1, function(x) identical(unname(x[1:n.var]), unname(dag.gprime[indices[[mcmc.move]][2],])))),ncol(sc.tmp)]
Q.G2 <- sumlogs(sc.tmp[,ncol(sc.tmp)])
##Orphaning
N.G <- sum(dag.tmp)
dag.gprime[indices[[mcmc.move]][1],] <- dag.gprime[indices[[mcmc.move]][2],] <- rep(0,n.var)
##step2:
sc.tmp <- sc[score.cache$children==indices[[mcmc.move]][2],]
sc.tmp <- sc.tmp[sc.tmp[,indices[[mcmc.move]][1]]==1,]
d.tmp <- dag.gprime
d.list.tmp <- NULL
s.list.tmp <-NULL
##test all b -> a
for(i in 1:dim(sc.tmp)[1]){
d.tmp[indices[[mcmc.move]][2],] <- sc.tmp[i,-ncol(sc.tmp)]
if(!identical(topoSortMAT(amat = d.tmp,index = FALSE),character(0))){
d.list.tmp[[length(d.list.tmp)+1]] <- d.tmp
s.list.tmp[[length(s.list.tmp)+1]] <- sc.tmp[i,ncol(sc.tmp)]
}}
step2 <- sample(x = 1:length(d.list.tmp),size = 1,replace = FALSE,prob = 1/abs(s.list.tmp))
##end step2:
dag.gprime[indices[[mcmc.move]][2],] <- sc.tmp[step2,-ncol(sc.tmp)]
Q.Gprime1 <- sumlogs(sc.tmp[,ncol(sc.tmp)])
##step3
sc.tmp <- sc[score.cache$children==indices[[mcmc.move]][1],]
d.tmp <- dag.gprime
d.list.tmp <- NULL
s.list.tmp <- NULL
##test all b -> a
for(i in 1:dim(sc.tmp)[1]){
d.tmp[indices[[mcmc.move]][1],] <- sc.tmp[i,-ncol(sc.tmp)]
if(!identical(topoSortMAT(amat = d.tmp,index = FALSE),character(0))){
d.list.tmp[[length(d.list.tmp)+1]] <- d.tmp
s.list.tmp[[length(s.list.tmp)+1]] <- sc.tmp[i,ncol(sc.tmp)]
}}
if(!is.null(s.list.tmp)){
step3 <- sample(x = 1:length(d.list.tmp),size = 1,replace = FALSE,prob = 1/abs(s.list.tmp))
##endstep3
dag.gprime[indices[[mcmc.move]][1],] <- sc.tmp[step3,-ncol(sc.tmp)]
Q.Gprime2 <- sumlogs(sc.tmp[,ncol(sc.tmp)])
N.Gprime <- sum(dag.gprime)
##output:
score.G <- score.Gprime <- 0
prior.G <- prior.Gprime <- 1
for(a in 1:n.var){
#Koivisto priors
if(prior==2){
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,]
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,])))),n.var+1]),score.Gprime)
}
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)]))
}
score <- score.G
#print(min(exp(log(prior.G) + log(N.G) + Q.G2 + Q.G1 - (log(prior.Gprime) + log(N.Gprime) + Q.Gprime2 + Q.Gprime1)),1))
#alpha<-1
alpha <- min(exp((log(prior.G) + log(N.G) + Q.G2 + Q.G1) - (log(prior.Gprime) + log(N.Gprime) + Q.Gprime2 + Q.Gprime1)),1)
}}#if tests!
##############################
## 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.j.M.cross)
A <- min(1,score.A)
if(rbinom(n = 1,size = 1,prob = A)==1){
if(rbinom(n = 1,size = 1,prob = alpha)==1){
dag.tmp <- dag.gprime
if(verbose)cat("REV move\n");
score <- score.Gprime
}
dag.tmp <- dag.M.tilde
score.M.tilde <- 0
for(a in 1:n.var){
return(list(dag.tmp=dag.tmp, score=score, alpha=alpha))
sc.tmp <- sc[score.cache$children==a,]
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
if(is.null(A)){A<-0}
##############################
## Return
##############################
return(list(dag.tmp = dag.tmp, score = score, alpha=A))
}#EOF
mc3<-function(n.var, dag.tmp, max.parents, sc,score.cache, score, prior, verbose){
mc3<-function(n.var, dag.tmp, max.parents, sc ,score.cache, score, prior, verbose){
##construction of neighbours list
neighbours.list <- NULL
......
......@@ -94,9 +94,9 @@ mcmc <- function(score.cache = NULL,
##REV move
####################################################
if(rbinom(n = 1,size = 1,prob = prob.rev)==1){
#if(verbose){print("REV move")}
out <- REV(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior,verbose)
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
......@@ -124,11 +124,9 @@ mcmc <- function(score.cache = NULL,
j <- 0
while(j<=mcmc.scheme[2]){
#if()
if(rbinom(n = 1,size = 1,prob = prob.rev)==1){
#if(verbose){print("fuck"); print( dag.tmp)}
out <- REV(n.var,(dag.tmp),max.parents,sc,score.cache,score,prior,verbose)
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
......
......@@ -126,3 +126,40 @@ formula.mcmcabn<-function(f, name){
##output
return(out)
}
##-------------------------------------------------------------------------
##Ancestor function
##-------------------------------------------------------------------------
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)
}
}
##-------------------------------------------------------------------------
##Descendant function
##-------------------------------------------------------------------------
descendants<-function(nodes, dag){
diag(dag)<-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{
tmp <- t(unname(which(x = dag[,nodes]==1,arr.ind = TRUE)))[1,]
}
return(unique(unname(c(tmp,descendants(nodes = (tmp),dag = dag)))))
}else{
return(NULL)
}
}
range01 <- function(x){(x-min(x))/(max(x)-min(x))}
......@@ -3,31 +3,31 @@ model{
###-----------------------
###Binomial nodes
###-----------------------
a ~ dbern(p.a); #Binary response
logit(p.a) <- 0.1 + 2.7*c; #Logistic regression
A ~ dbern(p.A); #Binary response
logit(p.A) <- 0 + 0.45*C; #Logistic regression
###-----------------------
###Binomial nodes
###-----------------------
b ~ dbern(p.b); #Binary response
logit(p.b) <- 0.1 + 2.7*a + 2.7*e; #Logistic regression
B ~ dbern(p.B); #Binary response
logit(p.B) <- 0.45 + 0.45*A + 0.45*E; #Logistic regression
###-----------------------
###Binomial nodes
###-----------------------
c ~ dbern(p.c); #Binary response
logit(p.c) <- 0.1; #Logistic regression
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 + 2.7*e; #Logistic regression
D ~ dbern(p.D); #Binary response
logit(p.D) <- 0 + 0.45*E; #Logistic regression
###-----------------------
###Binomial nodes
###-----------------------
e ~ dbern(p.e); #Binary response
logit(p.e) <- 0.1 + 2.7*c; #Logistic regression
E ~ dbern(p.E); #Binary response
logit(p.E) <- 0.45 + 0.45*C; #Logistic regression
}
\ No newline at end of file
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