Introduction

This document contains code to support the paper Margins of discrete Bayesian networks. The code finds the model defined by the observable margin of the Bayesian network model \(1 \rightarrow 2 \leftarrow U \rightarrow 4 \leftarrow 3\), when all variables (including the latent \(U\)) are binary.

Singular Code

We work with the computational algebra package Singular, which can be run in browser here.

We define parameters of the form pabcd to represent the probabilities \(P(X_2 = b, X_4 = d \mid X_1 = a, X_3 = c)\). These are written parametrically as \[ P(X_2 = b, X_4 = d \mid X_1 = a, X_3 = c) = \alpha q(b \mid a) q(d \mid c) + (1-\alpha) r(b \mid a) r(d \mid c), \] we then eliminate the parameters to obtain the ideal over the observed conditional probabilities.

LIB "sing.lib";

ring r = 0, (p0000,p1000,p0100,p1100,
p0010,p1010,p0110,p1110,
p0001,p1001,p0101,p1101,
p0011,p1011,p0111,p1111,
al, q2_0, q2_1, q4_0, q4_1,
r2_0, r2_1, r4_0, r4_1), dp;

ideal i = - p0000 + al*q2_0*q4_0 + (1-al)*r2_0*r4_0,
- p0100 + al*(1-q2_0)*q4_0 + (1-al)*(1-r2_0)*r4_0,
- p0001 + al*q2_0*(1-q4_0) + (1-al)*r2_0*(1-r4_0),
- p0101 + al*(1-q2_0)*(1-q4_0) + (1-al)*(1-r2_0)*(1-r4_0),
- p1000 + al*q2_1*q4_0 + (1-al)*r2_1*r4_0,
- p1100 + al*(1-q2_1)*q4_0 + (1-al)*(1-r2_1)*r4_0,
- p1001 + al*q2_1*(1-q4_0) + (1-al)*r2_1*(1-r4_0),
- p1101 + al*(1-q2_1)*(1-q4_0) + (1-al)*(1-r2_1)*(1-r4_0),
- p0010 + al*q2_0*q4_1 + (1-al)*r2_0*r4_1,
- p0110 + al*(1-q2_0)*q4_1 + (1-al)*(1-r2_0)*r4_1,
- p0011 + al*q2_0*(1-q4_1) + (1-al)*r2_0*(1-r4_1),
- p0111 + al*(1-q2_0)*(1-q4_1) + (1-al)*(1-r2_0)*(1-r4_1),
- p1010 + al*q2_1*q4_1 + (1-al)*r2_1*r4_1,
- p1110 + al*(1-q2_1)*q4_1 + (1-al)*(1-r2_1)*r4_1,
- p1011 + al*q2_1*(1-q4_1) + (1-al)*r2_1*(1-r4_1),
- p1111 + al*(1-q2_1)*(1-q4_1) + (1-al)*(1-r2_1)*(1-r4_1);

ideal j = eliminate(i, al*q2_0*q4_0*q2_1*q4_1*r2_0*r4_0*r2_1*r4_1);

j is the ideal defining the observable distributions, and it contains 9 polynomial constraints on the observed probabilities. These are:

j[1]=p0011-p1011+p0111-p1111
j[2]=p0001-p1001+p0101-p1101
j[3]=p1010+p1110+p1011+p1111-1
j[4]=p0010+p0110+p1011+p1111-1
j[5]=p1100-p1110+p1101-p1111
j[6]=p0100-p0110+p0101+p0011-p1011-p1111
j[7]=p1000+p1110+p1001+p1111-1
j[8]=p0000+p0110+p1001-p0101+p1101-p0011+p1011+p1111-1
j[9]=p1110*p0101*p1011-p0110*p1101*p1011-p1110*p1001*p0111-p1110*p1101*p0111-p1101*p1011*p0111+p0110*p1001*p1111+p1110*p0101*p1111+p0101*p1011*p1111-p1101*p0111*p1111+p0101*p1111^2+p1101*p0111-p0101*p1111

Since there are 16 conditional probabilities, this suggests that the model is 7 dimensional: we can check this by recreating the ideal in the ring without the additional parameters.

ring r2 = 0, (p0000,p1000,p0100,p1100,
p0010,p1010,p0110,p1110,
p0001,p1001,p0101,p1101,
p0011,p1011,p0111,p1111), dp;

/* use the constraints above */
ideal j2 = p0011-p1011+p0111-p1111,
p0001-p1001+p0101-p1101,
p1010+p1110+p1011+p1111-1,
p0010+p0110+p1011+p1111-1,
p1100-p1110+p1101-p1111,
p0100-p0110+p0101+p0011-p1011-p1111,
p1000+p1110+p1001+p1111-1,
p0000+p0110+p1001-p0101+p1101-p0011+p1011+p1111-1,
p1110*p0101*p1011-p0110*p1101*p1011-p1110*p1001*p0111-p1110*p1101*p0111-p1101*p1011*p0111+p0110*p1001*p1111+p1110*p0101*p1111+p0101*p1011*p1111-p1101*p0111*p1111+p0101*p1111^2+p1101*p0111-p0101*p1111;

/* this command finds the dimension */
dim(std(j2));

The answer is indeed 7.