// Luis David Garcia-Puente, last modified: 07.25.09 // Changsung Kang /////////////////////////////////////////////////////////////////////////////// version="$Id: Markov.lib,v 1.0.0.0 2009/07/25 18:09:30 levandov Exp$"; category="Algebraic Statistics"; info=" LIBRARY: Markov.lib Markov Relations for Bayesian Networks PROCEDURES: info(ideal) displays codim, degree, #mingens. bnet(int,intvec) creates a matrix. Bayesian net in topological order. nondec(int,intmat) nondescendents of a vertex pairMarkov(intmat) pairwise Markov property parent(int,intmat) parents of a vertex nondecminusparents(int,intmat) nondescendents(excluding the parents) of a vertex localMarkov(intmat) local Markov property subset(int,list) subset of a list children(int,intmat) children of a vertex Bayes_ball(list,list,intmat) computes a maximal set of vertices d-separated from a set given a conditioning set globalMarkov(intmat) global Markov property equivStatements(list,list) tests whether two conditional independence statements are equivalent next(intvec,int,intvec) index of the next variable in the ring sdec(intvec) decimal representation of a list probring(intvec,list) creates a probability distribution ring index(list,intvec) index of the corresponding indeterminate cartesian(list) cartesian product of elements in a list Pairs(list) set of all pairs of elements in a list levels(int) levels of a random variable Prob(list,intvec) marginalization over the subset of the random variables specified \"IND\" Quad(list,list,list,list,list,list,intvec) quadric associated to a fixed event StatementQuadrics(list,list,list,intvec) all quadrics associated to a conditional independence statement MarkovIdeal(list,intvec) computes the Markov ideal for a list of conditional independence statements torideal(ideal,intvec) computes the distinguished component of an ideal map_observable(def,intvec,int) creates a ring map from a ring associated to an observable probability distribution to basering "; LIB "general.lib"; LIB "elim.lib"; LIB "primdec.lib"; LIB "standard.lib"; LIB "presolve.lib"; /////////////////////////////////////////////////////////////////////////////// // The comments are still a little lousy but I'll work on it later. // General function to return information about the ideal proc info (ideal I) "USAGE: info(I); I ideal RETURN: list of integers a[1],a[2] and a[3] with: - a[1] the codimension of I - a[2] the degree of I - a[3] the number of minimal generators of I EXAMPLE: example info; shows an example" { int c = nvars(basering) - dim(I); int d = degree(I); int s = size(minbase(I)); list l = c,d,s; return(l); } example { "EXAMPLE:"; echo = 2; intvec d = 2,2,2,2; int n = size(d); def pdR = probring(d); setring pdR; intvec v = 1,1,0,0,1,1; intmat m = bnet(n,v); list l = localMarkov(m); ideal I = MarkovIdeal(l,d); info(I); } /* Definition of the matrix corresponding to the Bayesian Network. This procedure assumes that the DAG is in topological order n > n-1 > ... > 1 */ proc bnet (int n, intvec u) "USAGE: bnet(n,u); n int, u intvec RETURN: an n*n matrix whose lower triangle is given by u m[i,j] implies the existence of an edge (i,j) EXAMPLE: example bnet; shows an example" { int i,j,k; intmat m[n][n]; for (i=2; i<=n; i++) { for (j=1; j i) { l = insert(l,s,size(l)); } else { c = nd[j]; e = nondec(c,m); if (size(e) == 0) { l = insert(l,s,size(l)); } else { for (k=1; k<=size(e); k++) { if (e[k] == i) { check = 1; break; } } if (check == 0) { l = insert(l,s,size(l)); } else { e = delete(e,k); e2 = delete(nd,j); if (size(e) != size(e2)) { l = insert(l,s,size(l)); } else { if (size(e) != 0) { check = 0; for (k=1; k<=size(e); k++) { if (e[k] != e2[k]) { check = 1; break; } } if (check == 1) { l = insert(l,s,size(l)); } } } } } } } } } } return(l); } example { "EXAMPLE:"; echo = 2; intvec v = 1,1,0,0,1,1; intmat m = bnet(4,v); pairMarkov(m); } ////////////////////////////////////////////////////////////////// /* Definition of the Local Markov Relations */ ////////////////////////////////////////////////////////////////// proc parent (int v, intmat m) "USAGE: parent(v,m); n int, m intmat RETURN: list: the parents of the vertex v EXAMPLE: example parent; shows an example" { int n = ncols(m); list l; int i; for (i=1; i<=n; i++) { if (m[i,v] == 1) { l = insert(l,i,size(l)); } } return(l); } example { "EXAMPLE:"; echo = 2; intvec v = 1,1,0,0,1,1; intmat m = bnet(4,v); parent(1,m); } proc nondecminusparents (int v, intmat m) "USAGE: nondec(v,m); n int, m intmat RETURN: list: the nondescendents(excluding the parents) of the vertex v EXAMPLE: example nondecminusparents; shows an example" { int n = ncols(m); list l,s, visited; int i,j,k; s = v; for (i=1; i<=n; i++) { visited[i] = 0; } while (size(s) != 0) { k = s[1]; s = delete(s,1); visited[k] = 1; for (j=1; j<=n; j++) { if (m[k,j] == 1) { if (visited[j] == 0) { s = insert(s,j); } } } } for (i=1; i<=n; i++) { if (visited[i] == 0 and m[i,v] == 0) { l = insert(l,i,size(l)); } } return(l); } example { "EXAMPLE:"; echo = 2; intvec v = 1,1,0,0,1,1; intmat m = bnet(4,v); nondecminusparents(1,m); } proc localMarkov (intmat m) "USAGE: localMarkov(m); m intmat RETURN: l list: the local Markov property l[i] corresponds to the ith conditional independence statement, l[i][1] is independent of l[i][2] given l[i][3] EXAMPLE: example localMarkov; shows an example" { int n = ncols(m); list pa, s, l, nd, e; int i,c,check; for (i=1; i<=n; i++) { s = list(); s = insert(s,list(i)); nd = nondecminusparents(i,m); if (size(nd) != 0) { pa = parent(i,m); s = insert(s,nd,size(s)); s = insert(s,pa,size(s)); if (size(nd) > 1) { l = insert(l,s,size(l)); } else { c = nd[1]; if (c > i) { l = insert(l,s,size(l)); } else { e = parent(c,m); if (size(e) != size(pa)) { l = insert(l,s,size(l)); } else { if (size(e) != 0) { check = 1; while (check < size(e) and e[check] == pa[check]) { check++; } if (check != size(e)) { l = insert(l,s,size(l)); } } e = nondecminusparents(c,m); if (size(e) != 1) { l = insert(l,s,size(l)); } else { if (e[1] != i) { l = insert(l,s,size(l)); } } } } } } } return(l); } example { "EXAMPLE:"; echo = 2; intvec v = 1,1,0,0,1,1; intmat m = bnet(4,v); localMarkov(m); } ////////////////////////////////////////////////////////////////// /* Construction of the Global Markov Relations */ ////////////////////////////////////////////////////////////////// proc subset (int k, list X) "USAGE: subset(k,X); k int, X list RETURN: list: a subset of X @* If @math{b_n\cdots b_1} is the binary representation of the integer k, then subset(k,X) returns the set {@math{X[i]~|~b_i=1}} EXAMPLE: example subset; shows an example" { if (size(X) == 0) { return(list()); } int n = size(X); intvec bin; int q,i; for (i=1; i<=n; i++) { bin[i] = 0; } q = k; i = 1; while (q != 0) { bin[i] = q%2; q = q / 2; i++; } list res; for (i=1; i<=n; i++) { if (bin[i] == 1) { res = insert(res, X[i] ,size(res)); } } return(res); } example { "EXAMPLE:"; echo = 2; list l = 1,2,3; for(int i=1;i<=7;i++) { subset(i,l); } } proc children (int v, intmat m) "USAGE: children(v,m); n int, m intmat RETURN: list: the children of the vertex v EXAMPLE: example children; shows an example" { int n = ncols(m); list l; int i; for (i=1; i<=n; i++) { if (m[v,i] == 1) { l = insert(l,i,size(l)); } } return(l); } example { "EXAMPLE:"; echo = 2; intvec v = 1,1,0,0,1,1; intmat m = bnet(4,v); children(4,m); } proc Bayes_ball (list A, list C, intmat m) "USAGE: Bayes_ball(A,C,m); A list, C list, m intmat RETURN: list: a maximal set of vertices B such that A and B are d-separated by C EXAMPLE: example Bayes_ball; shows an example" { int n = ncols(m); int i,v; list B, pa,ch,vqueue; intvec visited, blocked, up, down, top, bottom; for (i=1; i<=n; i++) { visited[i] = 0; blocked[i] = 0; up[i] = 0; down[i] = 0; top[i] = 0; bottom[i] = 0; } for (i=1; i<=size(C); i++) { blocked[C[i]] = 1; } for (i=1; i<=size(A); i++) { vqueue = insert(vqueue, A[i]); up[A[i]] = 1; } while (size(vqueue) != 0) { v = vqueue[size(vqueue)]; vqueue = delete(vqueue,size(vqueue)); visited[v] = 1; if (!blocked[v] and up[v]) { if (!top[v]) { top[v] = 1; pa = parent(v,m); for (i=1; i<=size(pa); i++) { vqueue = insert(vqueue, pa[i]); up[pa[i]] = 1; } } if (!bottom[v]) { bottom[v] = 1; ch = children(v,m); for (i=1; i<=size(ch); i++) { vqueue = insert(vqueue, ch[i]); down[ch[i]] = 1; } } } if (down[v]) { if (blocked[v] and !top[v]) { top[v] = 1; pa = parent(v,m); for (i=1; i<=size(pa); i++) { vqueue = insert(vqueue, pa[i]); up[pa[i]] = 1; } } if (!blocked[v] and !bottom[v]) { bottom[v] = 1; ch = children(v,m); for (i=1; i<=size(ch); i++) { vqueue = insert(vqueue, ch[i]); down[ch[i]] = 1; } } } } for (i=1; i<=n; i++) { if(!bottom[i] and !blocked[i]) { B = insert(B,i,size(B)); } } return(B); } example { "EXAMPLE:"; echo = 2; intvec v = 1,0,1,0,1,0; intmat m = bnet(4,v); list A = 1; list C = 2; Bayes_ball(A,C,m); } proc globalMarkov (intmat m) "USAGE: globalMarkov(m); m intmat RETURN: l list: the global Markov property l[i] corresponds to the ith conditional independence statement, l[i][1] is independent of l[i][2] given l[i][3] EXAMPLE: example localMarkov; shows an example" { int n = ncols(m); int i,j,k,d,flag; list X,Y,A,B,C,l,s; for (i=1; i<=n; i++) { X[i] = i; } for (i=1; i<2^n-1; i++) { A = subset(i,X); Y = subset(2^n-i-1, X); d = size(Y); for (j=0; j<2^d-1; j++) { C = subset(j,Y); B = Bayes_ball(A,C,m); if (size(B)!=0) { flag = 0; s = list(); s = insert(s,A); s = insert(s,B,size(s)); s = insert(s,C,size(s)); for (k=1; k<=size(l); k++) { if (equivStatements(s,l[k])) { flag = 1; break; } } if (!flag) { l = insert(l,s,size(l)); } } } } return(l); } example { "EXAMPLE:"; echo = 2; intvec v = 1,1,0,0,1,1; intmat m = bnet(4,v); globalMarkov(m); } proc equivStatements(list s, list t) "USAGE: equivStatements(s,t); s list, t list RETURN: 1 if s[1]=t[2], s[2]=t[1] and s[3]=t[3] 0 otherwise EXAMPLE: example equivStatements; shows an example" { int i; if (size(s[1])!=size(t[2]) or size(s[2])!=size(t[1]) or size(s[3])!=size(t[3])) { return(0); } for (i=1; i<=size(s[1]); i++) { if (s[1][i] != t[2][i]) { return(0); } } for (i=1; i<=size(s[2]); i++) { if (s[2][i] != t[1][i]) { return(0); } } for (i=1; i<=size(s[3]); i++) { if (s[3][i] != t[3][i]) { return(0); } } return(1); } example { "EXAMPLE:"; echo = 2; list s = 1,list(2,3),4; list t = list(2,3),1,4; equivStatements(s,t); } ////////////////////////////////////////////////////////////////// /* Construction of the probability distribution ring */ ////////////////////////////////////////////////////////////////// /* This procedure computes the index of the next variable in the ring */ proc next (intvec u, int j, intvec d) "USAGE: next(u,j,d); u intvec, j int, d intvec RETURN: intvec: the index of the next variable in the ring EXAMPLE: example next; shows an example" { intvec v = u; if (j > 1) { v[j] = (v[j]+1)%(d[j]+1); if (v[j] == 0) { v[j] = 1; v = next (v, j-1, d); } } else { int check = (v[j]+1)%(d[j]+1); if (check != 0) { v[j] = check; } } return(v); } example { "EXAMPLE:"; echo = 2; intvec d = 2,2,2; def pdR = probring(d); setring pdR; intvec idx; for (int i=1; i<=size(d); i++) { idx[i] = 1; } idx = next(idx,size(d),d); idx; idx = next(idx,size(d),d); idx; } proc sdec (intvec id) "USAGE: sdec(id); id intvec RETURN: string: id[1]*10^(n-1)+id[2]*10^(n-2)+...+id[n] EXAMPLE: example sdec; shows an example" { int n = size(id); int dec; for (int i=n; i>=1; i--) { dec = dec + id[i]*(10)^(n-i); } return(string(dec)); } example { "EXAMPLE:"; echo = 2; intvec id = 1,4,10; sdec(id); } proc sdec2 (intvec id) { int n = size(id); string s; for (int i=1; i<=n; i++) { s = s+string(id[i]); } return(s); } /* This procedure has as its input the list of d_i's */ proc probring (intvec d, list#) "USAGE: probring(d[,f,v,o]); d intvec, f string, v string, o string RETURN: ring: ring R with coefficient field f, ring variables v1...1,...,vd[1]...d[n] and term ordering o @*The default values for f, v and o are \"0\", \"p\" and \"dp\" respectively EXAMPLE: example probring; shows an example" { if (size(#)==0 or size(#)>3) { #[1] = "0"; #[2] = "p"; #[3] = "dp";} if (size(#)==1) { #[2] = "p"; #[3] = "dp";} if (size(#)==2) { #[3] = "dp";} int i; intvec idx; for (i=1; i<=size(d); i++) { idx[i] = 1; } string ringconstructor = "ring nameofmyprettyring = " + #[1] + ",(" + #[2] + sdec(idx); for (i=2; i<=product(d); i++) { ringconstructor = ringconstructor + ","; idx = next(idx,size(d),d); ringconstructor = ringconstructor + #[2] + sdec(idx); } ringconstructor = ringconstructor + ")," + #[3] + ";"; execute(ringconstructor); return(basering); } example { "EXAMPLE:"; echo = 2; intvec d = 2,2,3; probring(d); } /* This procedure return the index of the corresponding indeterminate */ proc index (list linput, intvec d) "USAGE: index(linput,d); linput list, d intvec RETURN: int: the index of the corresponding indeterminate EXAMPLE: example index; shows an example" { int i; int base = 1; list shift = linput; for (i=1; i<=size(shift); i++) { shift[i] = shift[i]-1; } int n = size(shift); int idx = shift[n]; for (i=n-1; i>=1; i--) { base = base*d[i+1]; idx = idx + shift[i]*base; } idx = idx+1; return(idx); } example { "EXAMPLE:"; echo = 2; intvec d = 2,2,2; def pdR = probring(d); setring pdR; list l = 1,2,1; index(l,d); var(index(l,d)); } ////////////////////////////////////////////////////////////////// /* Computation of the ideal I_M for all the models defined above */ ////////////////////////////////////////////////////////////////// /* This procedure computes the cartesian product of a list of lists */ proc cartesian (list linput) "USAGE: cartesian(linput); linput list RETURN: list: the cartesian product of a list of lists EXAMPLE: example index; shows an example" { int i,j,k; if (size(linput) == 1) { list l = linput[1]; for (i=1; i<=size(l); i++) {l[i] = list(l[i]);} return(l); } list head = linput[1]; list tail = cartesian(delete(linput,1)); list final, each; for (i=1; i<=size(head); i++) { for (j=1; j<=size(tail); j++) { each = insert(tail[j], head[i]); final = insert(final, each, size(final)); } } return(final); } example { "EXAMPLE:"; echo = 2; list l = list(1,2),list(1,2),list(1); cartesian(l); } /* This procedure computes the set of all pairs of a given list */ proc Pairs (list L) "USAGE: Pairs(L); L list RETURN: list: the set of all pairs of L EXAMPLE: example Pairs; shows an example" { int i,j; list result; for (i=1; i<=size(L)-1; i++) { for (j=i+1; j<=size(L); j++) { result = insert(result, list(L[i],L[j]), size(result)); } } return(result); } example { "EXAMPLE:"; echo = 2; Pairs(list(1,2,3)); } /* This little procedure computes the levels of the random variable Xi */ proc levels (int di) "USAGE: levels(di); di int RETURN: list: the levels of the random variable Xi EXAMPLE: example levels; shows an example" { list l; for (int i=1; i<=di; i++) { l[i] = i; } return(l); } example { "EXAMPLE:"; echo = 2; levels(3); } /* This procedure computes the linear form Prob(A) for an instantiation of a subset A of the random variables */ proc Prob (list linput, intvec d) "USAGE: Prob(linput,d); linput list, d intvec RETURN: poly: the marginalization over the subset of the random variables specified \"IND\" EXAMPLE: example Prob; shows an example" { int i; list l=linput; for (i=1; i<=size(linput); i++) { if (typeof(l[i]) == "string") { l[i] = levels(d[i]); } else { l[i] = list(l[i]); } } l = cartesian(l); poly pr; for (i=1; i<=size(l); i++) { pr = pr + var(index(l[i],d)); } return(pr); } example { "EXAMPLE:"; echo = 2; intvec d = 2,2,2; def pdR = probring(d); setring pdR; list l = 1,"IND","IND"; Prob(l,d); } /* This procedure computes the quadric associated to the fixed event P(A=a,B=b,C=c) */ proc Quad (list A, list a, list B, list b, list C, list c, intvec d) "USAGE: Quad (A,a,B,b,C,c,d); A list, a list, B list, b list, C list, c list, d intvec RETURN: poly: the quadric associated to the probability P(A=a[1],B=b[1],C=c)* P(A=a[2],B=b[2],C=c)-P(A=a[2],B=b[1],C=c)*P(A=a[1],B=b[2],C=c) EXAMPLE: example Quad; shows an example" { int i,j; list Q1, Q2, Q3, Q4, l; poly P1,P2,P3,P4; for (i=1; i<=size(d); i++) { Q1[i] = "IND"; Q2[i] = "IND"; Q3[i] = "IND"; Q4[i] = "IND"; } if (size(C) != 0) { for (i=1; i<=size(C); i++) { j = C[i]; Q1[j]=c[i]; Q2[j]=c[i]; Q3[j]=c[i]; Q4[j]=c[i]; } } for (i=1; i<=size(B); i++) { j = B[i]; l = b[1]; Q1[j] = l[i]; Q3[j] = l[i]; l = b[2]; Q2[j] = l[i]; Q4[j] = l[i]; } for (i=1; i<=size(A); i++) { j = A[i]; l = a[1]; Q1[j] = l[i]; Q4[j] = l[i]; l = a[2]; Q2[j] = l[i]; Q3[j] = l[i]; } P1 = Prob(Q1,d); P2 = Prob(Q2,d); P3 = Prob(Q3,d); P4 = Prob(Q4,d); return(P1*P2 - P3*P4); } example { "EXAMPLE:"; echo = 2; /* Computes the probability P(X1=1,X2=1,X3=1)*P(X1=2,X2=2,X3=1) -P(X1=2,X2=1,X3=1)*P(X1=1,X2=2,X3=1) */ intvec d = 2,2,2; def pdR = probring(d); setring pdR; list A,B,C; list a,b,c; A = list(1); B = list(2); C = list(3); a[1] = levels(d[1]); b[1] = levels(d[2]); c[1] = levels(d[3]); a = Pairs(cartesian(a)); b = Pairs(cartesian(b)); c = cartesian(c); Quad(A,a[1],B,b[1],C,c[1],d); } /* This procedure computes the list of all quadrics associated to the probability of the event P(A,B,C) */ proc StatementQuadrics (list A, list B, list C, intvec d) "USAGE: StatementQuadrics(A,B,C,d); A list, B list, C list, d intvec RETURN: poly: the list of all quadrics associated to the conditional independence statement, A is independent of B given C EXAMPLE: example StatementQuadrics; shows an example" { int i,j,k; list a,b,c,result; for (i=1; i<=size(A); i++) { a[i] = levels(d[A[i]]); } a = Pairs(cartesian(a)); for (i=1; i<=size(B); i++) { b[i] = levels(d[B[i]]); } b = Pairs(cartesian(b)); if (size(C) == 0) { for (i=1; i<=size(a); i++) { for (j=1; j<=size(b); j++) { result = insert(result, Quad(A, a[i], B, b[j], C, list(), d), size(result)); } } } else { for (i=1; i<=size(C); i++) { c[i] = levels(d[C[i]]); } c = cartesian(c); for (k=1; k<=size(c); k++) { for (i=1; i<=size(a); i++) { for (j=1; j<=size(b); j++) { result = insert(result, Quad(A, a[i], B, b[j], C, c[k], d), size(result)); } } } } return(result); } example { "EXAMPLE:"; echo = 2; /* Lists all quadrics associated to the conditional independence statement, X1 is independent of X2 given X3 */ intvec d = 2,2,2; def pdR = probring(d); setring pdR; StatementQuadrics(list(1),list(2),list(3),d); } proc MarkovIdeal (list L, intvec d) "USAGE: MarkovIdeal(L,d); L list, d intvec RETURN: ideal: the ideal of the independence model given by L EXAMPLE: example MarkovIdeal; shows an example" { option(redSB); int i,j; ideal result; list l,genlist; for (i=1; i<=size(L); i++) { l = L[i]; genlist = StatementQuadrics(l[1],l[2],l[3],d); for (j=1; j<=size(genlist); j++) { result = result,genlist[j]; } result = interred(result); } // result = groebner(result); return(result); } example { "EXAMPLE:"; echo = 2; intvec d = 2,2,2,2; int n = size(d); def pdR = probring(d); setring pdR; intvec v15 = 1,1,0,0,1,1; intmat m15 = bnet(n,v15); list l15 = localMarkov(m15); list pw15 = pairMarkov(m15); list g15 = globalMarkov(m15); ideal I15 = MarkovIdeal(l15,d); info(I15); ideal G15 = MarkovIdeal(g15,d); info(G15); quotient(I15,G15); ideal T15 = torideal(I15,d); quotient(I15,T15); ideal Q15 = sat(I15,T15)[1]; list pd15 = primdecGTZ(Q15); info(T15)[1]; for (int i=1; i<=size(pd15); i++) { info(std(pd15[i][1]))[1]; } } proc torideal (ideal I, intvec d) "USAGE: torideal(I,d); I ideal, d intvec RETURN: ideal: I:p@math{^\infty} where p is the product of all the linear forms @*For example, if d=2,2,2, then p=p111*...p222*p+11*...p+22*p++1*p++2 EXAMPLE: example torideal; shows an example" { list l,f; int i,j; ideal t = I; for (i=1; i<=size(d); i++) { for (j=1; j 0) { p = quad[1]; pols = insert(pols,p); wv = findvars(lead(p)); m = m*wv[1]; mons = insert(mons,wv[2]); if (size(#) == 0) { J = eliminate(J,wv[1]); } else { J = fastelim(J,wv[1]); } quad = degreepart(shortid(J,2),2,2); } return(J,m,pols,mons); } proc fiber (ideal I, list pols, list mons) { ideal J = I; int i; for (i=1; i<=size(mons); i++) { J = sat(J+pols[i],mons[i])[1]; } return (J); } proc primdecGSS (ideal I, ideal T) { list components; ideal K = reduce(T,std(I)); K = K + 0; ideal Q; int c; for (c=1; c<=size(K); c++) { Q = sat(I,K[c])[1]; if (size(Q) > 1) { components = insert(components,Q); } } return (components); } proc hidden_preimage (ideal T, intvec d) { def t_ring = basering; intvec d2 = d[1..size(d)-1]; probring(d2, "h_ring", "p", "dp"); setring t_ring; int i,l; string s; poly form; ideal phi; for (i=1;i<=nvars(h_ring);i++) { s = "form = "; for (l=1; l<=d[size(d)]; l++) { s = s , varstr(h_ring,i) + string(l); if (l