Macaulay2 » Documentation
Packages » GraphicalModels :: gaussianParametrization
next | previous | forward | backward | up | index | toc

gaussianParametrization -- parametrization of the covariance matrix in terms of treks

Synopsis

Description

Given a mixed graph $G$ with directed and bidirected edges, let $L$ be the matrix corresponding to the directed edges (see directedEdgesMatrix) and let $W$ be the matrix corresponding to the bidirected edges (see bidirectedEdgesMatrix). Then, the covariance matrix $S$ (see covarianceMatrix) of the random variables in the Gaussian graphical model corresponding to the mixed graph $G$ can be parametrized by the matrix equation $S = (I-L)^{-T}W(I-L)^{-1}$, where $I$ is the identity matrix.

The entry $s_{(i,j)}$ of the covariance matrix can also be written as the sum of all monomials corresponding to treks between vertices $i$ and $j$. See trekSeparation for the definition of a trek. The monomial corresponding to a trek is the product of all parameters associated to the directed and bidirected edges on the trek.

The following example shows how to compute the ideal of the model using the parametrization, which could also be computed using gaussianVanishingIdeal

i1 : G = mixedGraph(digraph {{b,{c,d}},{c,{d}}},bigraph {{a,d}})

o1 = MixedGraph{Bigraph => Bigraph{a => {d}}   }
                                   d => {a}
                Digraph => Digraph{b => {c, d}}
                                   c => {d}
                                   d => {}
                Graph => Graph{}

o1 : MixedGraph
i2 : R = gaussianRing G

o2 = R

o2 : PolynomialRing
i3 : compactMatrixForm =false;
i4 : S = covarianceMatrix(R)

o4 = |  s     s     s     s     |
     |   a,a   a,b   a,c   a,d  |
     |                          |
     |  s     s     s     s     |
     |   a,b   b,b   b,c   b,d  |
     |                          |
     |  s     s     s     s     |
     |   a,c   b,c   c,c   c,d  |
     |                          |
     |  s     s     s     s     |
     |   a,d   b,d   c,d   d,d  |

             4      4
o4 : Matrix R  <-- R
i5 : L = directedEdgesMatrix(R)

o5 = |  0  0    0     0   |
     |                    |
     |  0  0  l     l     |
     |         b,c   b,d  |
     |                    |
     |  0  0    0   l     |
     |               c,d  |
     |                    |
     |  0  0    0     0   |

             4      4
o5 : Matrix R  <-- R
i6 : W = bidirectedEdgesMatrix(R)

o6 = |  p       0     0   p     |
     |   a,a               a,d  |
     |                          |
     |    0   p       0     0   |
     |         b,b              |
     |                          |
     |    0     0   p       0   |
     |               c,c        |
     |                          |
     |  p       0     0   p     |
     |   a,d               d,d  |

             4      4
o6 : Matrix R  <-- R
i7 : M = gaussianParametrization(R)

o7 = |  p                0                                0                  
     |   a,a                                                                 
     |                                                                       
     |    0             p                             l   p                  
     |                   b,b                           b,c b,b               
     |                                                                       
     |                                              2                        
     |    0           l   p                        l   p    + p              
     |                 b,c b,b                      b,c b,b    c,c           
     |                                                                       
     |                                  2                                    
     |  p     l   l   p    + l   p     l   l   p    + l   l   p    + l   p   
     |   a,d   b,c c,d b,b    b,d b,b   b,c c,d b,b    b,c b,d b,b    c,d c,c
     ------------------------------------------------------------------------
                                  p                                 |
                                   a,d                              |
                                                                    |
                        l   l   p    + l   p                        |
                         b,c c,d b,b    b,d b,b                     |
                                                                    |
                  2                                                 |
                 l   l   p    + l   l   p    + l   p                |
                  b,c c,d b,b    b,c b,d b,b    c,d c,c             |
                                                                    |
      2   2                              2          2               |
     l   l   p    + 2l   l   l   p    + l   p    + l   p    + p     |
      b,c c,d b,b     b,c b,d c,d b,b    b,d b,b    c,d c,c    d,d  |

             4      4
o7 : Matrix R  <-- R
i8 : J = delete(0_R, flatten entries (L|W))

o8 = {p   , p   , l   , l   , p   , l   , p   , p   , p   }
       a,a   a,d   b,c   b,d   b,b   c,d   c,c   a,d   d,d

o8 : List
i9 : eliminate(J, ideal(S-M))

o9 = ideal (s   , s   )
             a,c   a,b

o9 : Ideal of R
i10 : gaussianVanishingIdeal(R)

o10 = ideal (s   , s   )
              a,c   a,b

o10 : Ideal of R

This next example shows how to use the option SimpleTreks to compute a parametrization using simple treks instead of all treks. The resulting covariance matrix has diagonal entries equal to 1. This is giving a parametrization of all correlation matrices of matrices that belong to the model. This formulation is also known as Wright's method of path analysis.

i11 : G = mixedGraph(digraph {{b,{c,d}},{c,{d}}},bigraph {{a,d}})

o11 = MixedGraph{Bigraph => Bigraph{a => {d}}   }
                                    d => {a}
                 Digraph => Digraph{b => {c, d}}
                                    c => {d}
                                    d => {}
                 Graph => Graph{}

o11 : MixedGraph
i12 : R = gaussianRing G

o12 = R

o12 : PolynomialRing
i13 : M = gaussianParametrization(R,SimpleTreks=>true)

o13 = |    1          0                0               p          |
      |                                                 a,d       |
      |                                                           |
      |    0          1               l          l   l    + l     |
      |                                b,c        b,c c,d    b,d  |
      |                                                           |
      |    0         l                 1         l   l    + l     |
      |               b,c                         b,c b,d    c,d  |
      |                                                           |
      |  p     l   l    + l     l   l    + l            1         |
      |   a,d   b,c c,d    b,d   b,c b,d    c,d                   |

              4      4
o13 : Matrix R  <-- R

See also

Ways to use gaussianParametrization :

For the programmer

The object gaussianParametrization is a method function with options.