# objective: performing model determination for DAG's  with respect to a given ordering of the # variables, together with the corresponding parameter estimation

# required input: i) model specification; ii) data
# optional: iii) initial values


# i) model specification

# the complete model contains 5 binary vertices (A, B, C, D, E) and 10 directed arrows
# conditional independence relations are as follows
# pa(A)=empty
# pa(B)=A
# pa(C)=A, B
# pa(D)=A, B, C
# pa(E)=A, B, C, D

# declared variables:
# n1....                marginal counts for A
# n11...               marginal counts for B given A=1
# n01...               marginal counts for B given A=0
# nab1..              marginal counts for C given {A=a, B=b}
# nabc1.             marginal counts for D given {A=a, B=b, C=c}
# nabcd1            marginal counts for E given {A=a, B=b, C=c, D=d}

# the following appear in the paper as \epsilon_{A,B}, etc.
#eps[1]               indicator variable for arrow A->B (takes value 1 whenever the visited model 
#                          contains A->B ; 0 otherwise)
#eps[2]               indicator variable for arrow A->C
#eps[3]               indicator variable for arrow B->C
#eps[4]               indicator variable for arrow A->D
#eps[5]               indicator variable for arrow B->D
#eps[6]               indicator variable for arrow C->D
#eps[7]               indicator variable for arrow A->E
#eps[8]               indicator variable for arrow B->E
#eps[9]               indicator variable for arrow C->E
#eps[10]             indicator variable for arrow D->E

#the following four variables are used to compute the models posterior probabilities and will be explained below
# eta[j ]    ten-dimensional deterministic vector with components 2^j for j=9,...,0

# Z          unique number associated with each model 

# w[k]     index vector  for k=1,...,K (K=dimension of the model space, in this case 1024)

# I[k]        indicator variable for the kth model (takes value 1 whenever the corresponding         #              model is visited; 0 otherwise) 



#declared parameters:
# n                    number of trials for the A-marginal table
# n1....              number of trials for the B-marginal table given A=1
# n0....              number of trials for the B-marginal table given A=0
# nab...             number of trials for the C-marginal table given {A=a, B=b}
# nabc..            number of trials for the D-marginal table given {A=a, B=b, C=c}
# nabcd.           number of trials for the E-marginal table given {A=a, B=b, C=c, D=d}

# the following appear in the paper as \tilde{\theta^{V|pa(V)}_{v|pa(v)}}
# p1....              Pr{A=1}
# p11...             Pr{B=1|A=1}
# p01...             Pr{B=1|A=0}
# pab1..            Pr{C=1|A=a, B=b}
# pabc1.           Pr{D=1|A=a, B=b, C=c}
# pabcd1          Pr{E=1|A=a, B=b, C=c, D=d}
# for a generic vertex V,
# pVmarg         marginal probability of vertex V in {A, B, C, D, E} with the assumption that V      #                        has no parents
# pVWYwy        conditional probability of V given the combination of parents {W, Y} with 
#                        configuration {w, y}




model
{



#distribution of marginal counts: independent binomials, with cells probabilities expressed as #a linear combination of the conditional probabilities given each possible combination of #parents; for each realization of vector eps[ ] only one element of the combination is picked, #all the remaining ones being equal to zero

#vertex A:

n1....~dbin(p1....,n)
p1....<-pAmarg

#vertex B:

n11...~dbin(p11...,n1....)
p11...<-pBmarg*(1-eps[1])+pBA1*eps[1]

n0....<-n-n1....
n01...~dbin(p01...,n0....)
p01...<-pBmarg*(1-eps[1])+pBA0*eps[1]

#vertex C:

n111..~dbin(p111..,n11...)
p111..<- pCmarg*(1-eps[2])*(1-eps[3])+pCA1*eps[2]*(1-eps[3])+pCB1*eps[3]*(1-eps[2])+pCAB11*eps[2]*eps[3]

n10...<-n1....-n11...
n101..~dbin(p101..,n10...)
p101..<- pCmarg*(1-eps[2])*(1-eps[3])+pCA1*eps[2]*(1-eps[3])+pCB0*eps[3]*(1-eps[2])+pCAB10*eps[2]*eps[3]

n011..~dbin(p011..,n01...)
p011..<-pCmarg*(1-eps[2])*(1-eps[3])+pCA0*eps[2]*(1-eps[3])+pCB1*eps[3]*(1-eps[2])+pCAB01*eps[2]*eps[3]

n00...<-n0....-n01...
n001..~dbin(p001..,n00...)
p001..<-pCmarg*(1-eps[2])*(1-eps[3])+pCA0*eps[2]*(1-eps[3])+pCB0*eps[3]*(1-eps[2])+pCAB00*eps[2]*eps[3]

#vertex D:

n1111.~dbin(p1111.,n111..)
p1111.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA1*eps[4]*(1-eps[5])*(1-eps[6])+pDB1*(1-eps[4])*eps[5]*(1-eps[6])+pDC1*(1-eps[4])*(1-eps[5])*eps[6]+pDAB11*eps[4]*eps[5]*(1-eps[6])+pDAC11*eps[4]*(1-eps[5])*eps[6]+pDBC11*(1-eps[4])*eps[5]*eps[6]+pDABC111*eps[4]*eps[5]*eps[6]

n110..<-n11...-n111..
n1101.~dbin(p1101.,n110..)
p1101.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA1*eps[4]*(1-eps[5])*(1-eps[6])+pDB1*(1-eps[4])*eps[5]*(1-eps[6])+pDC0*(1-eps[4])*(1-eps[5])*eps[6]+pDAB11*eps[4]*eps[5]*(1-eps[6])+pDAC10*eps[4]*(1-eps[5])*eps[6]+pDBC10*(1-eps[4])*eps[5]*eps[6]+pDABC110*eps[4]*eps[5]*eps[6]

n1011.~dbin(p1011.,n101..)
p1011.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA1*eps[4]*(1-eps[5])*(1-eps[6])+pDB0*(1-eps[4])*eps[5]*(1-eps[6])+pDC1*(1-eps[4])*(1-eps[5])*eps[6]+pDAB10*eps[4]*eps[5]*(1-eps[6])+pDAC11*eps[4]*(1-eps[5])*eps[6]+pDBC01*(1-eps[4])*eps[5]*eps[6]+pDABC101*eps[4]*eps[5]*eps[6]

n0111.~dbin(p0111.,n011..)
p0111.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA0*eps[4]*(1-eps[5])*(1-eps[6])+pDB1*(1-eps[4])*eps[5]*(1-eps[6])+pDC1*(1-eps[4])*(1-eps[5])*eps[6]+pDAB01*eps[4]*eps[5]*(1-eps[6])+pDAC01*eps[4]*(1-eps[5])*eps[6]+pDBC11*(1-eps[4])*eps[5]*eps[6]+pDABC011*eps[4]*eps[5]*eps[6]

n0011.~dbin(p0011.,n001..)
p0011.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA0*eps[4]*(1-eps[5])*(1-eps[6])+pDB0*(1-eps[4])*eps[5]*(1-eps[6])+pDC1*(1-eps[4])*(1-eps[5])*eps[6]+pDAB00*eps[4]*eps[5]*(1-eps[6])+pDAC01*eps[4]*(1-eps[5])*eps[6]+pDBC01*(1-eps[4])*eps[5]*eps[6]+pDABC001*eps[4]*eps[5]*eps[6]

n010..<-n01...-n011..
n0101.~dbin(p0101.,n010..)
p0101.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA0*eps[4]*(1-eps[5])*(1-eps[6])+pDB1*(1-eps[4])*eps[5]*(1-eps[6])+pDC0*(1-eps[4])*(1-eps[5])*eps[6]+pDAB01*eps[4]*eps[5]*(1-eps[6])+pDAC00*eps[4]*(1-eps[5])*eps[6]+pDBC10*(1-eps[4])*eps[5]*eps[6]+pDABC010*eps[4]*eps[5]*eps[6]

n100..<-n10...-n101..
n1001.~dbin(p1001.,n100..)
p1001.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA1*eps[4]*(1-eps[5])*(1-eps[6])+pDB0*(1-eps[4])*eps[5]*(1-eps[6])+pDC0*(1-eps[4])*(1-eps[5])*eps[6]+pDAB10*eps[4]*eps[5]*(1-eps[6])+pDAC10*eps[4]*(1-eps[5])*eps[6]+pDBC00*(1-eps[4])*eps[5]*eps[6]+pDABC100*eps[4]*eps[5]*eps[6]

n000..<-n00...-n001..
n0001.~dbin(p0001.,n000..)
p0001.<-pDmarg*(1-eps[4])*(1-eps[5])*(1-eps[6])+pDA0*eps[4]*(1-eps[5])*(1-eps[6])+pDB0*(1-eps[4])*eps[5]*(1-eps[6])+pDC0*(1-eps[4])*(1-eps[5])*eps[6]+pDAB00*eps[4]*eps[5]*(1-eps[6])+pDAC00*eps[4]*(1-eps[5])*eps[6]+pDBC00*(1-eps[4])*eps[5]*eps[6]+pDABC000*eps[4]*eps[5]*eps[6]

#vertex E:

h89<-eps[9]*eps[10]
h.9<-(1-eps[9])*eps[10]
h8.<-eps[9]*(1-eps[10])
h..<-(1-eps[9])*(1-eps[10])


n11111~dbin(p11111,n1111.)

zero11111<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one11111<-pEA1*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two11111<- pEAB11*eps[7]*eps[8]*h..+pEAC11*eps[7]*(1-eps[8])*h8.+pEAD11*eps[7]*(1-eps[8])*h.9+pEBC11*(1-eps[7])*eps[8]*h8.+pEBD11*(1-eps[7])*eps[8]*h.9+pECD11*(1-eps[7])*(1-eps[8])*h89
three11111<-pEABC111*eps[7]*eps[8]*h8.+pEABD111*eps[7]*eps[8]*h.9+  pEACD111*eps[7]*(1-eps[8])*h89+pEBCD111*(1-eps[7])*eps[8]*h89  
four11111<-pEABCD1111*eps[7]*eps[8]*h89

p11111<-zero11111+one11111+two11111+three11111+four11111

n1110.<-n111..-n1111.
n11101~dbin(p11101,n1110.)

zero11101<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one11101<-pEA1*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two11101<- pEAB11*eps[7]*eps[8]*h..+pEAC11*eps[7]*(1-eps[8])*h8.+pEAD10*eps[7]*(1-eps[8])*h.9+pEBC11*(1-eps[7])*eps[8]*h8.+pEBD10*(1-eps[7])*eps[8]*h.9+pECD10*(1-eps[7])*(1-eps[8])*h89
three11101<-pEABC111*eps[7]*eps[8]*h8.+pEABD110*eps[7]*eps[8]*h.9+  pEACD110*eps[7]*(1-eps[8])*h89+pEBCD110*(1-eps[7])*eps[8]*h89  
four11101<-pEABCD1110*eps[7]*eps[8]*h89

p11101<-zero11101+one11101+two11101+three11101+four11101

n11011~dbin(p11011,n1101.)

zero11011<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one11011<-pEA1*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two11011<- pEAB11*eps[7]*eps[8]*h..+pEAC10*eps[7]*(1-eps[8])*h8.+pEAD11*eps[7]*(1-eps[8])*h.9+pEBC10*(1-eps[7])*eps[8]*h8.+pEBD11*(1-eps[7])*eps[8]*h.9+pECD01*(1-eps[7])*(1-eps[8])*h89
three11011<-pEABC110*eps[7]*eps[8]*h8.+pEABD111*eps[7]*eps[8]*h.9+  pEACD101*eps[7]*(1-eps[8])*h89+pEBCD101*(1-eps[7])*eps[8]*h89  
four11011<-pEABCD1101*eps[7]*eps[8]*h89

p11011<-zero11011+one11011+two11011+three11011+four11011

n10111~dbin(p10111,n1011.)

zero10111<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one10111<-pEA1*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two10111<- pEAB10*eps[7]*eps[8]*h..+pEAC11*eps[7]*(1-eps[8])*h8.+pEAD11*eps[7]*(1-eps[8])*h.9+pEBC01*(1-eps[7])*eps[8]*h8.+pEBD01*(1-eps[7])*eps[8]*h.9+pECD11*(1-eps[7])*(1-eps[8])*h89
three10111<-pEABC101*eps[7]*eps[8]*h8.+pEABD101*eps[7]*eps[8]*h.9+  pEACD111*eps[7]*(1-eps[8])*h89+pEBCD011*(1-eps[7])*eps[8]*h89  
four10111<-pEABCD1011*eps[7]*eps[8]*h89

p10111<-zero10111+one10111+two10111+three10111+four10111

n01111~dbin(p01111,n0111.)

zero01111<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one01111<-pEA0*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two01111<- pEAB01*eps[7]*eps[8]*h..+pEAC01*eps[7]*(1-eps[8])*h8.+pEAD01*eps[7]*(1-eps[8])*h.9+pEBC11*(1-eps[7])*eps[8]*h8.+pEBD11*(1-eps[7])*eps[8]*h.9+pECD11*(1-eps[7])*(1-eps[8])*h89
three01111<-pEABC011*eps[7]*eps[8]*h8.+pEABD011*eps[7]*eps[8]*h.9+  pEACD011*eps[7]*(1-eps[8])*h89+pEBCD111*(1-eps[7])*eps[8]*h89  
four01111<-pEABCD0111*eps[7]*eps[8]*h89

p01111<-zero01111+one01111+two01111+three01111+four01111

n00011~dbin(p00011,n0001.)

zero00011<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one00011<-pEA0*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two00011<- pEAB00*eps[7]*eps[8]*h..+pEAC00*eps[7]*(1-eps[8])*h8.+pEAD01*eps[7]*(1-eps[8])*h.9+pEBC00*(1-eps[7])*eps[8]*h8.+pEBD01*(1-eps[7])*eps[8]*h.9+pECD01*(1-eps[7])*(1-eps[8])*h89
three00011<-pEABC000*eps[7]*eps[8]*h8.+pEABD001*eps[7]*eps[8]*h.9+  pEACD001*eps[7]*(1-eps[8])*h89+pEBCD001*(1-eps[7])*eps[8]*h89  
four00011<-pEABCD0001*eps[7]*eps[8]*h89

p00011<-zero00011+one00011+two00011+three00011+four00011

n0010.<-n001..-n0011.
n00101~dbin(p00101,n0010.)

zero00101<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one00101<-pEA0*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two00101<- pEAB00*eps[7]*eps[8]*h..+pEAC01*eps[7]*(1-eps[8])*h8.+pEAD00*eps[7]*(1-eps[8])*h.9+pEBC01*(1-eps[7])*eps[8]*h8.+pEBD00*(1-eps[7])*eps[8]*h.9+pECD10*(1-eps[7])*(1-eps[8])*h89
three00101<-pEABC001*eps[7]*eps[8]*h8.+pEABD000*eps[7]*eps[8]*h.9+  pEACD010*eps[7]*(1-eps[8])*h89+pEBCD010*(1-eps[7])*eps[8]*h89  
four00101<-pEABCD0010*eps[7]*eps[8]*h89

p00101<-zero00101+one00101+two00101+three00101+four00101

n0100.<-n010..-n0101.
n01001~dbin(p01001,n0100.)

zero01001<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one01001<-pEA0*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two01001<- pEAB01*eps[7]*eps[8]*h..+pEAC00*eps[7]*(1-eps[8])*h8.+pEAD00*eps[7]*(1-eps[8])*h.9+pEBC10*(1-eps[7])*eps[8]*h8.+pEBD10*(1-eps[7])*eps[8]*h.9+pECD00*(1-eps[7])*(1-eps[8])*h89
three01001<-pEABC010*eps[7]*eps[8]*h8.+pEABD010*eps[7]*eps[8]*h.9+  pEACD000*eps[7]*(1-eps[8])*h89+pEBCD100*(1-eps[7])*eps[8]*h89  
four01001<-pEABCD0100*eps[7]*eps[8]*h89

p01001<-zero01001+one01001+two01001+three01001+four01001

n1000.<-n100..-n1001.
n10001~dbin(p10001,n1000.)

zero10001<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one10001<-pEA1*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two10001<- pEAB10*eps[7]*eps[8]*h..+pEAC10*eps[7]*(1-eps[8])*h8.+pEAD10*eps[7]*(1-eps[8])*h.9+pEBC00*(1-eps[7])*eps[8]*h8.+pEBD00*(1-eps[7])*eps[8]*h.9+pECD00*(1-eps[7])*(1-eps[8])*h89
three10001<-pEABC100*eps[7]*eps[8]*h8.+pEABD100*eps[7]*eps[8]*h.9+  pEACD100*eps[7]*(1-eps[8])*h89+pEBCD000*(1-eps[7])*eps[8]*h89  
four10001<-pEABCD1000*eps[7]*eps[8]*h89

p10001<-zero10001+one10001+two10001+three10001+four10001

n1010.<-n101..-n1011.
n10101~dbin(p10101,n1010.)

zero10101<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one10101<-pEA1*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two10101<- pEAB10*eps[7]*eps[8]*h..+pEAC11*eps[7]*(1-eps[8])*h8.+pEAD10*eps[7]*(1-eps[8])*h.9+pEBC01*(1-eps[7])*eps[8]*h8.+pEBD00*(1-eps[7])*eps[8]*h.9+pECD10*(1-eps[7])*(1-eps[8])*h89
three10101<-pEABC101*eps[7]*eps[8]*h8.+pEABD100*eps[7]*eps[8]*h.9+  pEACD110*eps[7]*(1-eps[8])*h89+pEBCD010*(1-eps[7])*eps[8]*h89  
four10101<-pEABCD1010*eps[7]*eps[8]*h89

p10101<-zero10101+one10101+two10101+three10101+four10101

n01011~dbin(p01011,n0101.)

zero01011<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one01011<-pEA0*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two01011<- pEAB01*eps[7]*eps[8]*h..+pEAC00*eps[7]*(1-eps[8])*h8.+pEAD01*eps[7]*(1-eps[8])*h.9+pEBC10*(1-eps[7])*eps[8]*h8.+pEBD11*(1-eps[7])*eps[8]*h.9+pECD01*(1-eps[7])*(1-eps[8])*h89
three01011<-pEABC010*eps[7]*eps[8]*h8.+pEABD011*eps[7]*eps[8]*h.9+  pEACD001*eps[7]*(1-eps[8])*h89+pEBCD101*(1-eps[7])*eps[8]*h89  
four01011<-pEABCD0101*eps[7]*eps[8]*h89

p01011<-zero01011+one01011+two01011+three01011+four01011

n10011~dbin(p10011,n1001.)

zero10011<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one10011<-pEA1*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two10011<- pEAB10*eps[7]*eps[8]*h..+pEAC10*eps[7]*(1-eps[8])*h8.+pEAD11*eps[7]*(1-eps[8])*h.9+pEBC00*(1-eps[7])*eps[8]*h8.+pEBD01*(1-eps[7])*eps[8]*h.9+pECD01*(1-eps[7])*(1-eps[8])*h89
three10011<-pEABC100*eps[7]*eps[8]*h8.+pEABD101*eps[7]*eps[8]*h.9+  pEACD101*eps[7]*(1-eps[8])*h89+pEBCD001*(1-eps[7])*eps[8]*h89  
four10011<-pEABCD1001*eps[7]*eps[8]*h89

p10011<-zero10011+one10011+two10011+three10011+four10011

n0110.<-n011..-n0111.
n01101~dbin(p01101,n0110.)

zero01101<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one01101<-pEA0*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two01101<- pEAB01*eps[7]*eps[8]*h..+pEAC01*eps[7]*(1-eps[8])*h8.+pEAD00*eps[7]*(1-eps[8])*h.9+pEBC11*(1-eps[7])*eps[8]*h8.+pEBD10*(1-eps[7])*eps[8]*h.9+pECD10*(1-eps[7])*(1-eps[8])*h89
three01101<-pEABC011*eps[7]*eps[8]*h8.+pEABD010*eps[7]*eps[8]*h.9+  pEACD010*eps[7]*(1-eps[8])*h89+pEBCD110*(1-eps[7])*eps[8]*h89  
four01101<-pEABCD0110*eps[7]*eps[8]*h89

p01101<-zero01101+one01101+two01101+three01101+four01101

n1100.<-n110..-n1101.
n11001~dbin(p11001,n1100.)

zero11001<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one11001<-pEA1*eps[7]*(1-eps[8])*h..+pEB1*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two11001<- pEAB11*eps[7]*eps[8]*h..+pEAC10*eps[7]*(1-eps[8])*h8.+pEAD10*eps[7]*(1-eps[8])*h.9+pEBC10*(1-eps[7])*eps[8]*h8.+pEBD10*(1-eps[7])*eps[8]*h.9+pECD00*(1-eps[7])*(1-eps[8])*h89
three11001<-pEABC110*eps[7]*eps[8]*h8.+pEABD110*eps[7]*eps[8]*h.9+  pEACD100*eps[7]*(1-eps[8])*h89+pEBCD100*(1-eps[7])*eps[8]*h89  
four11001<-pEABCD1100*eps[7]*eps[8]*h89

p11001<-zero11001+one11001+two11001+three11001+four11001

n00111~dbin(p00111,n0011.)

zero00111<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one00111<-pEA0*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC1*(1-eps[7])*(1-eps[8])*h8.+pED1*(1-eps[7])*(1-eps[8])*h.9
two00111<- pEAB00*eps[7]*eps[8]*h..+pEAC01*eps[7]*(1-eps[8])*h8.+pEAD01*eps[7]*(1-eps[8])*h.9+pEBC01*(1-eps[7])*eps[8]*h8.+pEBD01*(1-eps[7])*eps[8]*h.9+pECD11*(1-eps[7])*(1-eps[8])*h89
three00111<-pEABC001*eps[7]*eps[8]*h8.+pEABD001*eps[7]*eps[8]*h.9+  pEACD011*eps[7]*(1-eps[8])*h89+pEBCD011*(1-eps[7])*eps[8]*h89  
four00111<-pEABCD0011*eps[7]*eps[8]*h89

p00111<-zero00111+one00111+two00111+three00111+four00111

n0000.<-n000..-n0001.
n00001~dbin(p00001,n0000.)

zero00001<-pEmarg*(1-eps[7])*(1-eps[8])*h..
one00001<-pEA0*eps[7]*(1-eps[8])*h..+pEB0*(1-eps[7])*eps[8]*h..+pEC0*(1-eps[7])*(1-eps[8])*h8.+pED0*(1-eps[7])*(1-eps[8])*h.9
two00001<- pEAB00*eps[7]*eps[8]*h..+pEAC00*eps[7]*(1-eps[8])*h8.+pEAD00*eps[7]*(1-eps[8])*h.9+pEBC00*(1-eps[7])*eps[8]*h8.+pEBD00*(1-eps[7])*eps[8]*h.9+pECD00*(1-eps[7])*(1-eps[8])*h89
three00001<-pEABC000*eps[7]*eps[8]*h8.+pEABD000*eps[7]*eps[8]*h.9+  pEACD000*eps[7]*(1-eps[8])*h89+pEBCD000*(1-eps[7])*eps[8]*h89  
four00001<-pEABCD0000*eps[7]*eps[8]*h89

p00001<-zero00001+one00001+two00001+three00001+four00001


#prior distribution of conditional probability parameters: independent uniform

pAmarg ~ dunif(0,1)

 pBmarg~ dunif(0,1)
 pBA1~ dunif(0,1)
pBA0~ dunif(0,1)

   pCmarg~ dunif(0,1)
   pCB1~ dunif(0,1)
   pCB0~ dunif(0,1)
   pCA1~ dunif(0,1)
   pCA0~ dunif(0,1)
   pCAB11~ dunif(0,1)
   pCAB10~ dunif(0,1)
   pCAB01~ dunif(0,1)
   pCAB00~ dunif(0,1)

pDmarg~ dunif(0,1)
pDA1~ dunif(0,1)
pDA0~ dunif(0,1)
pDB1~ dunif(0,1)
pDB0~ dunif(0,1)
pDC1~ dunif(0,1)
pDC0~ dunif(0,1)
pDAB11~ dunif(0,1)
pDAB10~ dunif(0,1)
pDAB01~ dunif(0,1)
pDAB00~ dunif(0,1)
pDBC11~ dunif(0,1)
pDBC10~ dunif(0,1)
pDBC01~ dunif(0,1)
pDBC00~ dunif(0,1)
pDAC11~ dunif(0,1)
pDAC10~ dunif(0,1)
pDAC01~ dunif(0,1)
pDAC00~ dunif(0,1)
pDABC111~dunif(0,1)
pDABC110~dunif(0,1)
pDABC101~dunif(0,1)
pDABC011~dunif(0,1)
pDABC001~dunif(0,1)
pDABC010~dunif(0,1)
pDABC100~dunif(0,1)
pDABC000~dunif(0,1)

pEmarg~ dunif(0,1)
pEA1~ dunif(0,1)
pEA0~ dunif(0,1)
pEB1~ dunif(0,1)
pEB0~ dunif(0,1)
pEC1~ dunif(0,1)
pEC0~ dunif(0,1)
pED1~ dunif(0,1)
pED0~ dunif(0,1)
pEAB11~ dunif(0,1)
pEAB10~ dunif(0,1)
pEAB01~ dunif(0,1)
pEAB00~ dunif(0,1)
pEAC11~ dunif(0,1)
pEAC10~ dunif(0,1)
pEAC01~ dunif(0,1)
pEAC00~ dunif(0,1)
pEAD11~ dunif(0,1)
pEAD10~ dunif(0,1)
pEAD01~ dunif(0,1)
pEAD00~ dunif(0,1)
pEBC11~ dunif(0,1)
pEBC10~ dunif(0,1)
pEBC01~ dunif(0,1)
pEBC00~ dunif(0,1)
pEBD11~ dunif(0,1)
pEBD10~ dunif(0,1)
pEBD01~ dunif(0,1)
pEBD00~ dunif(0,1)
pECD11~ dunif(0,1)
pECD10~ dunif(0,1)
pECD01~ dunif(0,1)
pECD00~ dunif(0,1)
pEABC111~ dunif(0,1)
pEABC110~ dunif(0,1)
pEABC101~ dunif(0,1)
pEABC011~ dunif(0,1)
pEABC001~ dunif(0,1)
pEABC010~ dunif(0,1)
pEABC100~ dunif(0,1)
pEABC000~ dunif(0,1)
pEABD111~ dunif(0,1)
pEABD110~ dunif(0,1)
pEABD101~ dunif(0,1)
pEABD011~ dunif(0,1)
pEABD001~ dunif(0,1)
pEABD010~ dunif(0,1)
pEABD100~ dunif(0,1)
pEABD000~ dunif(0,1)
pEACD111~ dunif(0,1)
pEACD110~ dunif(0,1)
pEACD101~ dunif(0,1)
pEACD011~ dunif(0,1)
pEACD001~ dunif(0,1)
pEACD010~ dunif(0,1)
pEACD100~ dunif(0,1)
pEACD000~ dunif(0,1)
pEBCD111~ dunif(0,1)
pEBCD110~ dunif(0,1)
pEBCD101~ dunif(0,1)
pEBCD011~ dunif(0,1)
pEBCD001~ dunif(0,1)
pEBCD010~ dunif(0,1)
pEBCD100~ dunif(0,1)
pEBCD000~ dunif(0,1)
pEABCD1111~dunif(0,1)
pEABCD1110~dunif(0,1)
pEABCD1101~dunif(0,1)
pEABCD1011~dunif(0,1)
pEABCD0111~dunif(0,1)
pEABCD0001~dunif(0,1)
pEABCD0010~dunif(0,1)
pEABCD0100~dunif(0,1)
pEABCD1000~dunif(0,1)
pEABCD1010~dunif(0,1)
pEABCD0101~dunif(0,1)
pEABCD1001~dunif(0,1)
pEABCD0110~dunif(0,1)
pEABCD1100~dunif(0,1)
pEABCD0011~dunif(0,1)
pEABCD0000~dunif(0,1)




#prior for the arrow indicators: independent Bernoulli
   for( i in 1 : 10 ) {
      eps[i] ~ dbern(f[i])
   }


#indicator variable for each model
#notice that in order to calculate the posterior probabilities for each of the 1024 models we #have:
#a) identified each model with a unique number from 0 to 1023, depending of the realization #of vector eps [ ] and given by the result of the expansion inprod(e[],eta[ ]) (e.g. the complete model is identified by 1023; the model with no arrows by 0)
#b) associated to each of these numbers an indicator function I[k] which counts the number #of times each model appears in the entire MCMC run; by monitoring the components of I[ ] #one gets the posterior relative frequencies of all the models

   Z <- inprod(eps[],eta[ ])
   for( k in 1 : 1024 ) {
w[k]<-(k-1)
      I[k] <- equals(Z,w[k])
   }
}



# ii)data

list(eta=c(512,256,128,64,32,16,8,4,2,1),

	f=c(0.5,0.25,0.25,0.125,0.125,0.125,0.0625,0.0625,0.0625,0.0625),
	
		n=,
		
	n1....=,
	
	n11...=,
	n01...=,
	
	n111..=,
	n011..=,
	n101..=,
	n001..=,
	
	n1111.=,
	n1101.=,
	n1011.=,
	n0111.=,
	n0011.=,
	n0101.=,
	n1001.=,
	n0001.=,
	
	n11111=,
	n11101=,
	n11011=,
	n10111=,
	n01111=,
	n00011=,
	n00101=,
	n01001=,
	n10001=,
	n10101=,
	n01011=,
	n10011=,
	n01101=,
	n11001=,
	n00111=,
	n00001=

				
			
	)
	





# ii) initial values
   
list(eps=c(1,1,1,1,1,1,1,1,1,1),
pAmarg=0.5,
pBmarg=0.5,
pBA1=0.5,
pBA0=0.5,

pCmarg=0.5,
   pCA1=0.5,
   pCA0=0.5,
   pCB1=0.5,
   pCB0=0.5,
   pCAB11=0.5,
   pCAB10=0.5,
   pCAB01=0.5,
   pCAB00=0.5,

pDmarg=0.5,
pDB1=0.5,
pDB0=0.5,
pDC1=0.5,
pDC0=0.5,
pDA1=0.5,
pDA0=0.5,
pDBC11=0.5,
pDBC10=0.5,
pDBC01=0.5,
pDBC00=0.5,
pDAB11=0.5,
pDAB10=0.5,
pDAB01=0.5,
pDAB00=0.5,
pDAC11=0.5,
pDAC10=0.5,
pDAC01=0.5,
pDAC00=0.5,
pDABC111=0.5,
pDABC110=0.5,
pDABC101=0.5,
pDABC011=0.5,
pDABC001=0.5,
pDABC010=0.5,
pDABC100=0.5,
pDABC000=0.5,


pEmarg=0.5,
pEA1=0.5,
pEA0=0.5,
pEB1=0.5,
pEB0=0.5,
pEC1=0.5,
pEC0=0.5,
pED1=0.5,
pED0=0.5,
pEAB11=0.5,
pEAB10=0.5,
pEAB01=0.5,
pEAB00=0.5,
pEAC11=0.5,
pEAC10=0.5,
pEAC01=0.5,
pEAC00=0.5,
pEAD11=0.5,
pEAD10=0.5,
pEAD01=0.5,
pEAD00=0.5,
pEBC11=0.5,
pEBC10=0.5,
pEBC01=0.5,
pEBC00=0.5,
pEBD11=0.5,
pEBD10=0.5,
pEBD01=0.5,
pEBD00=0.5,
pECD11=0.5,
pECD10=0.5,
pECD01=0.5,
pECD00=0.5,
pEABC111=0.5,
pEABC110=0.5,
pEABC101=0.5,
pEABC011=0.5,
pEABC001=0.5,
pEABC010=0.5,
pEABC100=0.5,
pEABC000=0.5,
pEABD111=0.5,
pEABD110=0.5,
pEABD101=0.5,
pEABD011=0.5,
pEABD001=0.5,
pEABD010=0.5,
pEABD100=0.5,
pEABD000=0.5,

pEACD111=0.5,
pEACD110=0.5,
pEACD101=0.5,
pEACD011=0.5,
pEACD001=0.5,
pEACD010=0.5,
pEACD100=0.5,
pEACD000=0.5,
pEBCD111=0.5,
pEBCD110=0.5,
pEBCD101=0.5,
pEBCD011=0.5,
pEBCD001=0.5,
pEBCD010=0.5,
pEBCD100=0.5,
pEBCD000=0.5,
pEABCD1111=0.5,
pEABCD1110=0.5,
pEABCD1101=0.5,
pEABCD1011=0.5,
pEABCD0111=0.5,
pEABCD0001=0.5,
pEABCD0010=0.5,
pEABCD0100=0.5,
pEABCD1000=0.5,
pEABCD1010=0.5,
pEABCD0101=0.5,
pEABCD1001=0.5,
pEABCD0110=0.5,
pEABCD1100=0.5,
pEABCD0011=0.5,
pEABCD0000=0.5)