-- -*- coding: utf-8-unix -*-
-*
Copyright 2020 Carlos Amendola, Luis David Garcia Puente, Roser Homs Pons,
Olga Kuznetsova, Harshit J Motwani, Elina Robeva and David Swinarski.
You may redistribute this file under the terms of the GNU General Public
License as published by the Free Software Foundation, either version 2 of
the License, or any later version.
*-
newPackage(
"GraphicalModelsMLE",
Version => "1.0",
Date => "April 19, 2022",
Authors => {
{Name=> "Carlos Amendola",
Email=> "carlos.amendola@mis.mpg.de",
HomePage=>"http://www.carlos-amendola.com/"},
{Name => "Luis David Garcia Puente",
Email => "lgarciapuente@coloradocollege.edu",
HomePage => "https://luisgarciapuente.github.io"},
{Name=> "Roser Homs Pons",
Email=> "roser.homs@tum.de",
HomePage=>"https://personal-homepages.mis.mpg.de/homspons/index.html"},
{Name=> "Olga Kuznetsova",
Email=> "kuznetsova.olga@gmail.com",
HomePage=>"https://okuznetsova.com"},
{Name=> "Harshit J Motwani",
Email=> "harshitmotwani2015@gmail.com",
HomePage=>"https://sites.google.com/view/harshitjmotwani/home"},
{Name=> "Elina Robeva",
Email=> "erobeva@gmail.com",
HomePage=>"https://www.math.ubc.ca/~erobeva/"},
{Name=> "David Swinarski",
Email=> "dswinarski@fordham.edu",
HomePage=>"http://faculty.fordham.edu/dswinarski"}
},
Headline => "maximum likelihood estimates for graphical statistical models",
Keywords => {"Algebraic Statistics"},
PackageExports => {"GraphicalModels","Graphs","EigenSolver","NumericalAlgebraicGeometry","StatGraphs"},
Certification => {
"journal name" => "The Journal of Software for Algebra and Geometry",
"journal URI" => "https://msp.org/jsag/",
"article title" => "Computing maximum likelihood estimates for Gaussian graphical models with Macaulay2",
"acceptance date" => "14 March 2022",
"published article URI" => "https://msp.org/jsag/2022/12-1/p01.xhtml",
"published article DOI" => "10.2140/jsag.2022.12.1",
"published code URI" => "https://msp.org/jsag/2022/12-1/jsag-v12-n1-x01-GraphicalModelsMLE.zip",
"repository code URI" => "https://github.com/Macaulay2/M2/blob/master/M2/Macaulay2/packages/GraphicalModelsMLE.m2",
"release at publication" => "f0aaf5e74216283f022f7773bc1116f8de6b944b", -- git commit number in hex
"version at publication" => "1.0",
"volume number" => "12",
"volume URI" => "https://msp.org/jsag/2022/12-1/"
}
)
export {
"checkPD",
"checkPSD",
"Solver",--optional argument in solverMLE
"ConcentrationMatrix",-- optional argument in solverMLE
"CovarianceMatrix", -- optional argument in scoreEquations
"Saturate",-- optional argument in scoreEquations and solverMLE
"jacobianMatrixOfRationalFunction",
"MLdegree",
"OptionsEigenSolver",--optional argument in solverMLE
"OptionsNAG4M2",--optional argument in solverMLE
"RealPrecision",--optional argument in solverMLE and scoreEquations
"sampleCovarianceMatrix",
"SampleData",-- optional argument in scoreEquations and solverMLE
"SaturateOptions", -- optional argument in scoreEquations and solverMLE
"scoreEquations",
"solverMLE",
"ZeroTolerance"
}
--**************************--
-- INTERNAL ROUTINES --
--**************************--
--*************************************--
-- Functions (local) used throughout --
--*************************************--
------------------------------------------------------
-- Substitutes a list of points on a list of matrices
-- input - list of points from sols
-- matrix whose entries are variables
-- (expect it to be an inverse of a covariance matrix, Sin)
-- output - list of matrices after substituting these values
------------------------------------------------------
genListMatrix = (L,A) ->
(
T:= for l in L list coordinates(l);
M:= for t in T list substitute(A,matrix{t});
return M
);
----------------------------------------------
-- Selects all argmax for log det K- trace(S*K),
-- where K is an element of the list L.
-- We assume that L is the intersection of the
-- variety of the ideal generated by the Jacobian
-- of critical equations and the cone of PD matrices.
-- input - list L of candidate Sinv matrices (Sinv is Sigma^{-1}) and
-- sample covariance matrix V. Notation in line with scoreEquationsFromCovarianceMatrix
-- output -list of argmax
----------------------------------------------
maxMLE=(L,V)->(
if #L==0 then error("No critical points to evaluate");
if #L==1 then (E:=L_0; maxPt:=log det L_0- trace (V*L_0))
else
(eval:=for Sinv in L list log det Sinv- trace (V*Sinv);
evalReal:=for pt in eval when isReal pt list pt;
if #evalReal==0 then error("No critical point evaluates to a real solution");
maxPt=max evalReal;
indexOptimal:=positions(eval, i ->i== maxPt);
E= for i in indexOptimal list L_i;);
return (maxPt, E)
);
-------------------------------------------
-- scoreEquationsInternal - function that returns
-- both the ideal and the corresponding SInv matrix.
-- The user-facing scoreEquations method returns only
-- the ideal, whereas SInv is used in solverMLE
-------------------------------------------
scoreEquationsInternal={Saturate => true, SaturateOptions => options saturate, SampleData=>true, RealPrecision=>53, CovarianceMatrix => false}>>opts->(R,U)->(
----------------------------------------------------
-- Extract information about the graph
----------------------------------------------------
-- Lambda
L := directedEdgesMatrix R;
-- Psi
P := bidirectedEdgesMatrix R;
-- If the mixedGraph only has undirected part, call specific function for undirected.
if L==0 and P==0 then
return scoreEquationsInternalUndir(R,U,opts);
-- K
K := undirectedEdgesMatrix R;
----------------------------------------------------
-- Create an auxiliary ring and its fraction field
-- which do not have the s variables
----------------------------------------------------
-- create a new ring, lpR, which does not have the s variables
lpR:=coefficientRing(R)[gens R-set support covarianceMatrix R];
-- create its fraction field
FR := frac(lpR);
-----------------------------------------------------
-- Compute Sinv
-----------------------------------------------------
-- Kinv
K=sub(K, FR);
Kinv:=inverse K;
P=sub(P,FR);
--Omega
if K==0 then W:=P else (if P==0 then W=Kinv else W = directSum(Kinv,P));
-- move to FR, the fraction field of lpR
L= sub(L,FR);
-- Sigma
d:=numcols L;
if L==0 then S:=W else (
IdL := inverse (id_(FR^d)-L);
S = (transpose IdL) * W * IdL
);
Sinv := inverse S;
-----------------------------------------------------
-- Compute score equations ideal
----------------------------------------------------
-- Sample covariance matrix
if opts.SampleData then V:= sampleCovarianceMatrix(U) else V=U;
if ring V===RR_53 then V = matrix apply(entries V,r->r/(v->lift(numeric(opts.RealPrecision,v),QQ)));
-- Jacobian of log-likelihood function
C1 := trace(Sinv * V);
C1derivative := jacobianMatrixOfRationalFunction(C1);
LL :=jacobianMatrixOfRationalFunction (det Sinv)*matrix{{1/det(Sinv)}} - C1derivative;
LL=flatten entries(LL);
denoms := apply(#LL, i -> lift(denominator(LL_i), lpR));
J:=ideal apply(#LL, i -> lift(numerator(LL_i),lpR));
--Saturate
if opts.Saturate then (
argSaturate:=opts.SaturateOptions >>newOpts-> args ->(args, newOpts);
for i from 0 to (#denoms-1) do (
if degree denoms_i =={0} then J=J else
J=saturate(argSaturate(J,denoms_i))
);
);
return (J,Sinv);
);
----------------------------------------------------
--scoreEquationsInternalUndir for undirected graphs
----------------------------------------------------
scoreEquationsInternalUndir={Saturate => true, SaturateOptions => options saturate, SampleData=>true, RealPrecision=> 53, CovarianceMatrix => false}>>opts->(R,U)->(
-- Sample covariance matrix
if opts.SampleData then V := sampleCovarianceMatrix(U) else V=U;
if ring V===RR_53 then V = matrix apply(entries V,r->r/(v->lift(numeric(opts.RealPrecision,v),QQ)));
-- Concentration matrix K
K:=undirectedEdgesMatrix R;
-- move to a new ring, lpR, which does not have the s variables
lpR:=coefficientRing(R)[gens R - set support covarianceMatrix R];
K=sub(K,lpR);
J:=ideal{jacobian ideal{determinant(K)}-determinant(K)*jacobian(ideal{trace(K*V)})};
if opts.Saturate then
( argSaturate:=opts.SaturateOptions >>newOpts-> args ->(args, newOpts);
J=saturate(argSaturate(J,ideal{determinant(K)}));
);
return (J,K);
);
-------------------------------------------
-- Method copied from package DeterminantalRepresentations
-------------------------------------------
-----------------------------------------------
-- Method for retrieving the real part of a matrix.
--The code of this function is directly taken from DeterminantalRepresentations package in M2.
realPartMatrix = A -> matrix apply(entries A, r -> r/realPart)
--**************************--
-- METHODS --
--**************************--
sampleCovarianceMatrix = method(TypicalValue =>Matrix);
sampleCovarianceMatrix(Matrix) := (U) -> (
n := numRows U;
--Convert from integers to rationals if needed
if ring U===ZZ then U=sub(U,QQ);
--Convert matrix into list of row matrices
U = for i to n-1 list U^{i};
--Compute the mean vector
Ubar := matrix{{(1/n)}} * sum(U);
--Compute sample covariance matrix
return ((1/n)*(sum apply(n, i -> (transpose (U#i-Ubar))*(U#i-Ubar))));
);
sampleCovarianceMatrix(List) := (U) -> (
return sampleCovarianceMatrix(matrix U);
);
jacobianMatrixOfRationalFunction = method(TypicalValue =>Matrix);
jacobianMatrixOfRationalFunction(RingElement) := (F) -> (
if not instance(ring F,FractionField) then error "Expected element in a field of fractions";
f:=numerator(F);
g:=denominator(F);
R:=ring(f);
answer:=diff(vars(R), f) * g - diff(vars(R), g)*f;
answer=substitute(answer, ring(F));
return transpose(matrix({{(1/g)^2}})*answer)
);
scoreEquations = method(TypicalValue =>Sequence, Options =>{SampleData => true, Saturate => true, SaturateOptions => options saturate, RealPrecision => 53, CovarianceMatrix => false});
scoreEquations(Ring,Matrix) := opts -> (R, U) -> (
----------------------------------------------------
--Check input
----------------------------------------------------
if not R.?graph then error "Expected a ring created with gaussianRing of a Graph, Bigraph, Digraph or MixedGraph";
if not numRows U==#vertices R.graph then error "Size of sample data does not match the graph.";
if not opts.SampleData then (if not U==transpose U then error "The sample covariance matrix must be symmetric.");
---------------------------------------------------
-- Apply appropriate scoreEquations routine
---------------------------------------------------
if R.graphType===Graph
then (J,Sinv):=scoreEquationsInternalUndir(R,U,opts)
else (J,Sinv)=scoreEquationsInternal(R,U,opts);
if not opts.CovarianceMatrix
then return J
else return (J, inverse Sinv)
;
);
scoreEquations(Matrix,Ring) := opts ->(U,R) -> (
return scoreEquations(R,U,opts);
);
scoreEquations(Ring,List) := opts ->(R, U) -> (
----------------------------------------------------
--Check input
----------------------------------------------------
if not opts.SampleData then error "The sample covariance matrix must be a matrix.";
---------------------------------------------------
-- Call scoreEquations routine with a matrix
---------------------------------------------------
return scoreEquations(R,matrix U,opts);
);
scoreEquations(List,Ring) := opts ->(U,R) -> (
return scoreEquations(R,U,opts);
);
--Allow for graphs as input
scoreEquations(MixedGraph, Matrix) := opts ->(G,U) -> (
return scoreEquations(gaussianRing(G),U,opts);
);
scoreEquations(MixedGraph,List) := opts ->(G,U) -> (
return scoreEquations(gaussianRing(G),U,opts);
);
--All permutations of input for MixedGraphs and other type of graphs
scoreEquations(Graph,List) := opts -> (G, U) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Digraph,List) := opts -> (G, U) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Bigraph,List) := opts -> (G, U) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Graph,Digraph,List) := opts -> (G,D,U) -> (
return scoreEquations(mixedGraph (G,D),U, opts);
);
scoreEquations(Digraph,Graph,List) := opts -> (D,G,U) -> (
return scoreEquations(mixedGraph (D,G),U, opts);
);
scoreEquations(Digraph,Bigraph,List) := opts -> (D,B,U) -> (
return scoreEquations(mixedGraph (D,B),U, opts);
);
scoreEquations(Bigraph,Digraph,List) := opts -> (B,D,U) -> (
return scoreEquations(mixedGraph (B,D),U, opts);
);
scoreEquations(Graph, Bigraph,List) := opts -> (G,B,U) -> (
return scoreEquations(mixedGraph (G,B),U, opts);
);
scoreEquations(Bigraph,Graph,List) := opts -> (B,G,U) -> (
return scoreEquations(mixedGraph (B,G),U, opts);
);
scoreEquations(Graph, Digraph, Bigraph, List) := opts -> (G,D,B,U) -> (
return scoreEquations(mixedGraph (G,D,B),U, opts);
);
scoreEquations(Digraph, Bigraph, Graph, List) := opts -> (D,B,G,U) -> (
return scoreEquations(mixedGraph (D,B,G),U, opts);
);
scoreEquations(Bigraph, Graph, Digraph, List) := opts -> (B,G,D,U) -> (
return scoreEquations(mixedGraph (B,G,D),U, opts);
);
scoreEquations(Graph,Bigraph, Digraph, List) := opts -> (G,B,D,U) -> (
return scoreEquations(mixedGraph (G,B,D),U, opts);
);
scoreEquations(Bigraph, Digraph,Graph, List) := opts -> (B,D,G,U) -> (
return scoreEquations(mixedGraph (B,D,G),U, opts);
);
scoreEquations(Digraph, Graph, Bigraph, List) := opts -> (D,G,B,U) -> (
return scoreEquations(mixedGraph (D,G,B),U, opts);
);
scoreEquations(List,MixedGraph) := opts -> (U,G) -> (
return scoreEquations(G,U, opts);
);
scoreEquations(List,Graph) := opts -> (U,G) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(List,Digraph) := opts -> (U,G) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(List,Bigraph) := opts -> (U,G) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(List,Graph,Digraph) := opts -> (U,G,D) -> (
return scoreEquations(mixedGraph (G,D),U, opts);
);
scoreEquations(List,Digraph,Graph) := opts -> (U,D,G) -> (
return scoreEquations(mixedGraph (D,G),U, opts);
);
scoreEquations(List,Digraph,Bigraph) := opts -> (U,D,B) -> (
return scoreEquations(mixedGraph (D,B),U, opts);
);
scoreEquations(List,Bigraph,Digraph) := opts -> (U,B,D) -> (
return scoreEquations(mixedGraph (B,D),U, opts);
);
scoreEquations(List,Graph, Bigraph) := opts -> (U,G,B) -> (
return scoreEquations(mixedGraph (G,B),U, opts);
);
scoreEquations(List,Bigraph,Graph) := opts -> (U,B,G) -> (
return scoreEquations(mixedGraph (B,G),U, opts);
);
scoreEquations(List,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
return scoreEquations(mixedGraph (G,D,B),U, opts);
);
scoreEquations(List,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
return scoreEquations(mixedGraph (D,B,G),U, opts);
);
scoreEquations(List,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
return scoreEquations(mixedGraph (B,G,D),U, opts);
);
scoreEquations(List,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
return scoreEquations(mixedGraph (G,B,D),U, opts);
);
scoreEquations(List,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
return scoreEquations(mixedGraph (B,D,G),U, opts);
);
scoreEquations(List,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
return scoreEquations(mixedGraph (D,G,B),U, opts);
);
scoreEquations(Graph,Matrix) := opts -> (G, U) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Digraph,Matrix) := opts -> (G, U) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Bigraph,Matrix) := opts -> (G, U) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Graph,Digraph,Matrix) := opts -> (G,D,U) -> (
return scoreEquations(mixedGraph (G,D),U, opts);
);
scoreEquations(Digraph,Graph,Matrix) := opts -> (D,G,U) -> (
return scoreEquations(mixedGraph (D,G),U, opts);
);
scoreEquations(Digraph,Bigraph,Matrix) := opts -> (D,B,U) -> (
return scoreEquations(mixedGraph (D,B),U, opts);
);
scoreEquations(Bigraph,Digraph,Matrix) := opts -> (B,D,U) -> (
return scoreEquations(mixedGraph (B,D),U, opts);
);
scoreEquations(Graph, Bigraph,Matrix) := opts -> (G,B,U) -> (
return scoreEquations(mixedGraph (G,B),U, opts);
);
scoreEquations(Bigraph,Graph,Matrix) := opts -> (B,G,U) -> (
return scoreEquations(mixedGraph (B,G),U, opts);
);
scoreEquations(Graph, Digraph, Bigraph, Matrix) := opts -> (G,D,B,U) -> (
return scoreEquations(mixedGraph (G,D,B),U, opts);
);
scoreEquations(Digraph, Bigraph, Graph, Matrix) := opts -> (D,B,G,U) -> (
return scoreEquations(mixedGraph (D,B,G),U, opts);
);
scoreEquations(Bigraph, Graph, Digraph, Matrix) := opts -> (B,G,D,U) -> (
return scoreEquations(mixedGraph (B,G,D),U, opts);
);
scoreEquations(Graph,Bigraph, Digraph, Matrix) := opts -> (G,B,D,U) -> (
return scoreEquations(mixedGraph (G,B,D),U, opts);
);
scoreEquations(Bigraph, Digraph,Graph, Matrix) := opts -> (B,D,G,U) -> (
return scoreEquations(mixedGraph (B,D,G),U, opts);
);
scoreEquations(Digraph, Graph, Bigraph, Matrix) := opts -> (D,G,B,U) -> (
return scoreEquations(mixedGraph (D,G,B),U, opts);
);
scoreEquations(Matrix,MixedGraph) := opts -> (U,G) -> (
return scoreEquations(G,U, opts);
);
scoreEquations(Matrix,Graph) := opts -> (U,G) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Matrix,Digraph) := opts -> (U,G) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Matrix,Bigraph) := opts -> (U,G) -> (
return scoreEquations(mixedGraph (G),U, opts);
);
scoreEquations(Matrix,Graph,Digraph) := opts -> (U,G,D) -> (
return scoreEquations(mixedGraph (G,D),U, opts);
);
scoreEquations(Matrix,Digraph,Graph) := opts -> (U,D,G) -> (
return scoreEquations(mixedGraph (D,G),U, opts);
);
scoreEquations(Matrix,Digraph,Bigraph) := opts -> (U,D,B) -> (
return scoreEquations(mixedGraph (D,B),U, opts);
);
scoreEquations(Matrix,Bigraph,Digraph) := opts -> (U,B,D) -> (
return scoreEquations(mixedGraph (B,D),U, opts);
);
scoreEquations(Matrix,Graph, Bigraph) := opts -> (U,G,B) -> (
return scoreEquations(mixedGraph (G,B),U, opts);
);
scoreEquations(Matrix,Bigraph,Graph) := opts -> (U,B,G) -> (
return scoreEquations(mixedGraph (B,G),U, opts);
);
scoreEquations(Matrix,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
return scoreEquations(mixedGraph (G,D,B),U, opts);
);
scoreEquations(Matrix,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
return scoreEquations(mixedGraph (D,B,G),U, opts);
);
scoreEquations(Matrix,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
return scoreEquations(mixedGraph (B,G,D),U, opts);
);
scoreEquations(Matrix,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
return scoreEquations(mixedGraph (G,B,D),U, opts);
);
scoreEquations(Matrix,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
return scoreEquations(mixedGraph (B,D,G),U, opts);
);
scoreEquations(Matrix,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
return scoreEquations(mixedGraph (D,G,B),U, opts);
);
checkPD = method(TypicalValue =>List, Options =>{ZeroTolerance=>1e-10});
checkPD(List) := opts -> (L) -> (
for l in L
list (
if not length (select(eigenvalues l, i->realPart i<= opts.ZeroTolerance
or abs(imaginaryPart i )>opts.ZeroTolerance))==0 then continue;
realPartMatrix l)
);
checkPD(Matrix):= opts -> (L)->{
return checkPD({L},opts);
};
checkPSD = method(TypicalValue =>List, Options =>{ZeroTolerance=>1e-10});
checkPSD(List) := opts -> (L) -> (
for l in L
list (
if not length (select(eigenvalues l, i->realPart i< -opts.ZeroTolerance
or abs(imaginaryPart i )>opts.ZeroTolerance))==0 then continue;
realPartMatrix l)
);
checkPSD(Matrix):= opts -> (L)->{
return checkPSD({L},opts);
};
MLdegree = method(TypicalValue =>ZZ);
MLdegree(Ring):= (R) -> (
if not R.?graph then error "Expected gaussianRing created from a graph, digraph, bigraph or mixedGraph";
n:=# vertices R.graph;
J:=scoreEquations(R,random(RR^n,RR^n));
dimJ := dim J;
if dimJ > 0 then error concatenate("the ideal of score equations has dimension ",toString dimJ, " > 0,
so ML degree is not well-defined. The degree of this ideal is ", toString degree J,".");
return degree J;
);
solverMLE = method(TypicalValue =>Sequence, Options =>{SampleData=>true, ConcentrationMatrix=> false, Saturate => true, SaturateOptions => options saturate, Solver=>"EigenSolver", OptionsEigenSolver => options zeroDimSolve, OptionsNAG4M2=> options solveSystem, RealPrecision => 53, ZeroTolerance=>1e-10});
solverMLE(MixedGraph,Matrix) := opts -> (G, U) -> (
return solverMLE(gaussianRing G,U,opts);
);
-- Allow list instead of matrix
solverMLE(MixedGraph,List):= opts ->(G,U) -> (
return solverMLE(gaussianRing G,U,opts);
);
--Allow ring instead of graph
solverMLE(Ring,Matrix):= opts ->(R,U) -> (
-- check input
if not numgens source U ==R.gaussianRingData#nn then error "Size of sample data does not match the graph.";
-- sample covariance matrix
if opts.SampleData then V := sampleCovarianceMatrix(U)
else (V=U;
if not V==transpose V then error "The sample covariance matrix must be symmetric.");
-- generate the ideal of the score equations
if opts.Saturate then (
argSaturate:=opts.SaturateOptions >>newOpts-> args ->(args, SaturateOptions=>newOpts,SampleData=>false);
(J,SInv):=scoreEquationsInternal(argSaturate(R,V));)
else (J,SInv)= scoreEquationsInternal(R,V,Saturate=>false, SampleData=>false);
-- check that the system has finitely many solutions
if dim J =!= 0 then return J
else (
ML:=degree J;
-- solve system
if opts.Solver=="EigenSolver" then(
argES:=opts.OptionsEigenSolver >>newOpts-> args ->(args, newOpts);
sols:=zeroDimSolve(argES(J));
) else (
if opts.Solver=="NAG4M2" then (
sys:= flatten entries gens J;
argNAG4M2:=opts.OptionsNAG4M2 >>newOpts-> args ->(args, newOpts);
sols=solveSystem(argNAG4M2(sys));
)
else error "Accepted solver options are EigenSolver
(which uses function zeroDimSolve) or NAG4M2 (which uses solveSystem). Options should
be given as strings.";
);
--evaluate matrices on solutions
M:=genListMatrix(sols,SInv);
--consider only PD matrices
L:=checkPD (M, ZeroTolerance=>opts.ZeroTolerance);
--find the optimal points
(maxPt, E):=maxMLE(L,V);
if not opts.ConcentrationMatrix then (
if instance(E,List) then E=(for e in E list e=inverse e) else E=inverse E
);
return (maxPt,E,ML));
);
-- Allow ring instead of graph and list instead of matrix
solverMLE(Ring,List):= opts ->(R,U) -> (
-- check input
if not opts.SampleData then error "The sample covariance matrix must be a matrix.";
-- call solverMLE for a matrix
return solverMLE(R,matrix U,opts);
);
-- Permutations of input
solverMLE(Matrix,Ring) := opts -> (U, R) -> (
return solverMLE(R,U, opts);
);
solverMLE(List,Ring) := opts -> (U, R) -> (
return solverMLE(R,U, opts);
);
solverMLE(Graph,List) := opts -> (G, U) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Digraph,List) := opts -> (G, U) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Bigraph,List) := opts -> (G, U) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Graph,Digraph,List) := opts -> (G,D,U) -> (
return solverMLE(mixedGraph (G,D),U, opts);
);
solverMLE(Digraph,Graph,List) := opts -> (D,G,U) -> (
return solverMLE(mixedGraph (D,G),U, opts);
);
solverMLE(Digraph,Bigraph,List) := opts -> (D,B,U) -> (
return solverMLE(mixedGraph (D,B),U, opts);
);
solverMLE(Bigraph,Digraph,List) := opts -> (B,D,U) -> (
return solverMLE(mixedGraph (B,D),U, opts);
);
solverMLE(Graph, Bigraph,List) := opts -> (G,B,U) -> (
return solverMLE(mixedGraph (G,B),U, opts);
);
solverMLE(Bigraph,Graph,List) := opts -> (B,G,U) -> (
return solverMLE(mixedGraph (B,G),U, opts);
);
solverMLE(Graph, Digraph, Bigraph, List) := opts -> (G,D,B,U) -> (
return solverMLE(mixedGraph (G,D,B),U, opts);
);
solverMLE(Digraph, Bigraph, Graph, List) := opts -> (D,B,G,U) -> (
return solverMLE(mixedGraph (D,B,G),U, opts);
);
solverMLE(Bigraph, Graph, Digraph, List) := opts -> (B,G,D,U) -> (
return solverMLE(mixedGraph (B,G,D),U, opts);
);
solverMLE(Graph,Bigraph, Digraph, List) := opts -> (G,B,D,U) -> (
return solverMLE(mixedGraph (G,B,D),U, opts);
);
solverMLE(Bigraph, Digraph,Graph, List) := opts -> (B,D,G,U) -> (
return solverMLE(mixedGraph (B,D,G),U, opts);
);
solverMLE(Digraph, Graph, Bigraph, List) := opts -> (D,G,B,U) -> (
return solverMLE(mixedGraph (D,G,B),U, opts);
);
solverMLE(List,MixedGraph) := opts -> (U,G) -> (
return solverMLE(G,U, opts);
);
solverMLE(List,Graph) := opts -> (U,G) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(List,Digraph) := opts -> (U,G) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(List,Bigraph) := opts -> (U,G) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(List,Graph,Digraph) := opts -> (U,G,D) -> (
return solverMLE(mixedGraph (G,D),U, opts);
);
solverMLE(List,Digraph,Graph) := opts -> (U,D,G) -> (
return solverMLE(mixedGraph (D,G),U, opts);
);
solverMLE(List,Digraph,Bigraph) := opts -> (U,D,B) -> (
return solverMLE(mixedGraph (D,B),U, opts);
);
solverMLE(List,Bigraph,Digraph) := opts -> (U,B,D) -> (
return solverMLE(mixedGraph (B,D),U, opts);
);
solverMLE(List,Graph, Bigraph) := opts -> (U,G,B) -> (
return solverMLE(mixedGraph (G,B),U, opts);
);
solverMLE(List,Bigraph,Graph) := opts -> (U,B,G) -> (
return solverMLE(mixedGraph (B,G),U, opts);
);
solverMLE(List,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
return solverMLE(mixedGraph (G,D,B),U, opts);
);
solverMLE(List,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
return solverMLE(mixedGraph (D,B,G),U, opts);
);
solverMLE(List,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
return solverMLE(mixedGraph (B,G,D),U, opts);
);
solverMLE(List,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
return solverMLE(mixedGraph (G,B,D),U, opts);
);
solverMLE(List,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
return solverMLE(mixedGraph (B,D,G),U, opts);
);
solverMLE(List,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
return solverMLE(mixedGraph (D,G,B),U, opts);
);
solverMLE(Graph,Matrix) := opts -> (G, U) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Digraph,Matrix) := opts -> (G, U) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Bigraph,Matrix) := opts -> (G, U) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Graph,Digraph,Matrix) := opts -> (G,D,U) -> (
return solverMLE(mixedGraph (G,D),U, opts);
);
solverMLE(Digraph,Graph,Matrix) := opts -> (D,G,U) -> (
return solverMLE(mixedGraph (D,G),U, opts);
);
solverMLE(Digraph,Bigraph,Matrix) := opts -> (D,B,U) -> (
return solverMLE(mixedGraph (D,B),U, opts);
);
solverMLE(Bigraph,Digraph,Matrix) := opts -> (B,D,U) -> (
return solverMLE(mixedGraph (B,D),U, opts);
);
solverMLE(Graph, Bigraph,Matrix) := opts -> (G,B,U) -> (
return solverMLE(mixedGraph (G,B),U, opts);
);
solverMLE(Bigraph,Graph,Matrix) := opts -> (B,G,U) -> (
return solverMLE(mixedGraph (B,G),U, opts);
);
solverMLE(Graph, Digraph, Bigraph, Matrix) := opts -> (G,D,B,U) -> (
return solverMLE(mixedGraph (G,D,B),U, opts);
);
solverMLE(Digraph, Bigraph, Graph, Matrix) := opts -> (D,B,G,U) -> (
return solverMLE(mixedGraph (D,B,G),U, opts);
);
solverMLE(Bigraph, Graph, Digraph, Matrix) := opts -> (B,G,D,U) -> (
return solverMLE(mixedGraph (B,G,D),U, opts);
);
solverMLE(Graph,Bigraph, Digraph, Matrix) := opts -> (G,B,D,U) -> (
return solverMLE(mixedGraph (G,B,D),U, opts);
);
solverMLE(Bigraph, Digraph,Graph, Matrix) := opts -> (B,D,G,U) -> (
return solverMLE(mixedGraph (B,D,G),U, opts);
);
solverMLE(Digraph, Graph, Bigraph, Matrix) := opts -> (D,G,B,U) -> (
return solverMLE(mixedGraph (D,G,B),U, opts);
);
solverMLE(Matrix,MixedGraph) := opts -> (U,G) -> (
return solverMLE(G,U, opts);
);
solverMLE(Matrix,Graph) := opts -> (U,G) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Matrix,Digraph) := opts -> (U,G) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Matrix,Bigraph) := opts -> (U,G) -> (
return solverMLE(mixedGraph (G),U, opts);
);
solverMLE(Matrix,Graph,Digraph) := opts -> (U,G,D) -> (
return solverMLE(mixedGraph (G,D),U, opts);
);
solverMLE(Matrix,Digraph,Graph) := opts -> (U,D,G) -> (
return solverMLE(mixedGraph (D,G),U, opts);
);
solverMLE(Matrix,Digraph,Bigraph) := opts -> (U,D,B) -> (
return solverMLE(mixedGraph (D,B),U, opts);
);
solverMLE(Matrix,Bigraph,Digraph) := opts -> (U,B,D) -> (
return solverMLE(mixedGraph (B,D),U, opts);
);
solverMLE(Matrix,Graph, Bigraph) := opts -> (U,G,B) -> (
return solverMLE(mixedGraph (G,B),U, opts);
);
solverMLE(Matrix,Bigraph,Graph) := opts -> (U,B,G) -> (
return solverMLE(mixedGraph (B,G),U, opts);
);
solverMLE(Matrix,Graph, Digraph, Bigraph) := opts -> (U,G,D,B) -> (
return solverMLE(mixedGraph (G,D,B),U, opts);
);
solverMLE(Matrix,Digraph, Bigraph, Graph) := opts -> (U,D,B,G) -> (
return solverMLE(mixedGraph (D,B,G),U, opts);
);
solverMLE(Matrix,Bigraph, Graph, Digraph) := opts -> (U,B,G,D) -> (
return solverMLE(mixedGraph (B,G,D),U, opts);
);
solverMLE(Matrix,Graph,Bigraph, Digraph) := opts -> (U,G,B,D) -> (
return solverMLE(mixedGraph (G,B,D),U, opts);
);
solverMLE(Matrix,Bigraph, Digraph,Graph) := opts -> (U,B,D,G) -> (
return solverMLE(mixedGraph (B,D,G),U, opts);
);
solverMLE(Matrix,Digraph, Graph, Bigraph) := opts -> (U,D,G,B) -> (
return solverMLE(mixedGraph (D,G,B),U, opts);
);
--******************************************--
-- DOCUMENTATION --
--******************************************--
beginDocumentation()
doc ///
Key
GraphicalModelsMLE
Headline
a package for MLE of parameters for Gaussian graphical models
Description
Text
{\bf Graphical Models MLE} is a package for algebraic statistics that broadens the functionalities of @TO GraphicalModels@.
It computes the maximum likelihood estimates (MLE) of the covariance matrix of Gaussian graphical models associated to loopless mixed graphs(LMG).
The main features of the package are the computation of the @TO sampleCovarianceMatrix@ of sample data,
the ideal generated by @TO scoreEquations@ of log-likelihood functions of Gaussian graphical model,
the @TO MLdegree@ of such models and the MLE for the covariance or concentration matrix via @TO solverMLE@.
For more details on the type of graphical models that are accepted see @TO gaussianRing@.
In particular, for further information about LMG with undirected, directed and bidirected edges, check @TO partitionLMG@.
{\bf References:}
An introduction to key notions such as MLE and ML-degree can be found in the books:
Seth Sullivant, {\em Algebraic statistics}, American Mathematical Society, Vol 194, 2018.
Mathias Drton, Bernd Sturmfels and Seth Sullivant, {\em Lectures on Algebraic Statistics}, Oberwolfach Seminars, Vol 40, Birkhauser, Basel, 2009.
The definition and classification of loopless mixed graphs (LMG) can be found in the paper:
Kayvan Sadeghi and Steffen Lauritzen, {\em Markov properties for mixed graphs}, Bernoulli, 20 (2014), no 2, 676-696.
{\bf Examples:}
Computation of a sample covariance matrix from sample data:
Example
U= matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
sampleCovarianceMatrix(U)
Text
The ideal generated by the score equations of the log-likelihood function of the graphical model associated to the
graph $1\rightarrow 2,1\rightarrow 3,2\rightarrow 3,3\rightarrow 4,3<-> 4$ is computed as follows:
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
R = gaussianRing(G);
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
scoreEquations(R,U)
Text
Computation of the ML-degree of the 4-cycle:
Example
G=graph{{1,2},{2,3},{3,4},{4,1}};
MLdegree(gaussianRing G)
Text
Next compute the MLE for the covariance matrix of the graphical model associated
to the graph $1\rightarrow 3,2\rightarrow 4,3<-> 4,1 - 2$.
The input is the sample covariance instead of the sample data.
Example
G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}},graph {{1,2}});
V = matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}};
solverMLE(G,V,SampleData=>false)
Text
As an application of @TO solverMLE@: positive definite matrix completion
Text
Consider the following symmetric matrix with some unknown entries:
Example
R=QQ[x,y];
M=matrix{{115,-13,x,47},{-13,5,7,y},{x,7,27,-21},{47,y,-21,29}}
Text
Unknown entries correspond to the non-edges of the 4-cycle. A positive definite completion of this matrix
is obtained by giving values to x and y and computing the MLE for the covariance matrix in the Gaussian graphical model
given by the 4-cycle. Check @TO solverMLE@ for more details.
Example
G=graph{{1,2},{2,3},{3,4},{1,4}};
V=matrix{{115,-13,-29,47},{-13,5,7,-11},{-29,7,27,-21},{47,-11,-21,29}};
(mx,MLE,ML)=solverMLE(G,V,SampleData=>false)
Caveat
GraphicalModelsMLE requires @TO Graphs@, @TO StatGraphs@ and @TO GraphicalModels@.
In order to use the default numerical solver, it also requires @TO EigenSolver@.
@TO Graphs@ allows the user to create graphs whose vertices are labeled arbitrarily.
However, several functions in GraphicalModels sort the vertices of the graph. Hence, graphs used as input to methods
in GraphicalModelsMLE must have sortable vertex labels, e.g., all numbers or all letters.
@TO StatGraphs@ allows the user to work with objects such as bigraphs and mixedGraphs.
@TO GraphicalModels@ is used to generate @TO gaussianRing@, i.e. rings encoding graph properties.
///
--------------------------------
-- Documentation
--------------------------------
doc ///
Key
sampleCovarianceMatrix
(sampleCovarianceMatrix, List)
(sampleCovarianceMatrix, Matrix)
Headline
sample covariance matrix of observation vectors
Usage
sampleCovarianceMatrix U
Inputs
U:Matrix
or @TO List@ of sample data
Outputs
:Matrix
sample covariance matrix of the sample data
Description
Text
The sample covariance matrix is $S = \frac{1}{n} \sum_{i=1}^{n} (X^{(i)}-\bar{X}) (X^{(i)}-\bar{X})^T$.
Note that for normally distributed random variables, $S$ is the maximum likelihood estimator (MLE) for the
covariance matrix. This is different from the unbiased estimator, which uses a denominator of $n-1$ instead of $n$.
Sample data is input as a matrix or a list.
The rows of the matrix or the elements of the list are observation vectors.
Example
L= {{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
sampleCovarianceMatrix(L)
U= matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
sampleCovarianceMatrix(U)
///
doc ///
Key
jacobianMatrixOfRationalFunction
(jacobianMatrixOfRationalFunction,RingElement)
Headline
Jacobian matrix of a rational function
Usage
jacobianMatrixOfRationalFunction(F)
Inputs
F:RingElement
in @TO frac@
Outputs
:Matrix
the Jacobian matrix of a rational function
Description
Text
This function computes the Jacobian matrix of a rational function.
The input is an element in a fraction field.
Example
R=QQ[x,y];
FR=frac R;
F=1/(x^2+y^2);
jacobianMatrixOfRationalFunction(F)
Example
R=QQ[t_1,t_2,t_3];
FR=frac R;
jacobianMatrixOfRationalFunction( (t_1^2*t_2)/(t_1+t_2^2+t_3^3) )
///
-------------------------------------------------------
-- Documentation scoreEquations -----------------------
-------------------------------------------------------
doc ///
Key
scoreEquations
(scoreEquations, Ring, List)
(scoreEquations, Ring, Matrix)
(scoreEquations, List, Ring)
(scoreEquations, Matrix, Ring)
(scoreEquations, MixedGraph, Matrix)
(scoreEquations, MixedGraph, List)
(scoreEquations, Graph, List)
(scoreEquations, Digraph, List)
(scoreEquations, Bigraph, List)
(scoreEquations, Graph, Digraph,List)
(scoreEquations, Digraph, Graph,List)
(scoreEquations, Digraph, Bigraph, List)
(scoreEquations, Bigraph, Digraph, List)
(scoreEquations, Graph, Bigraph, List)
(scoreEquations, Bigraph, Graph, List)
(scoreEquations, Graph, Digraph, Bigraph, List)
(scoreEquations, Digraph, Bigraph, Graph, List)
(scoreEquations, Bigraph, Graph, Digraph, List)
(scoreEquations, Graph, Bigraph, Digraph, List)
(scoreEquations, Bigraph, Digraph, Graph, List)
(scoreEquations, Digraph, Graph, Bigraph, List)
(scoreEquations, List, MixedGraph)
(scoreEquations, List, Graph)
(scoreEquations, List, Digraph)
(scoreEquations, List, Bigraph)
(scoreEquations, List, Graph, Digraph)
(scoreEquations, List, Digraph,Graph)
(scoreEquations, List, Digraph,Bigraph)
(scoreEquations, List, Bigraph,Digraph)
(scoreEquations, List, Graph, Bigraph)
(scoreEquations, List, Bigraph, Graph)
(scoreEquations, List, Graph, Digraph, Bigraph)
(scoreEquations, List, Digraph, Bigraph, Graph)
(scoreEquations, List, Bigraph, Graph, Digraph)
(scoreEquations, List, Graph, Bigraph, Digraph)
(scoreEquations, List, Bigraph, Digraph, Graph)
(scoreEquations, List, Digraph, Graph, Bigraph)
(scoreEquations, Graph, Matrix)
(scoreEquations, Digraph, Matrix)
(scoreEquations, Bigraph, Matrix)
(scoreEquations, Graph, Digraph,Matrix)
(scoreEquations, Digraph, Graph,Matrix)
(scoreEquations, Digraph, Bigraph, Matrix)
(scoreEquations, Bigraph, Digraph, Matrix)
(scoreEquations, Graph, Bigraph, Matrix)
(scoreEquations, Bigraph, Graph, Matrix)
(scoreEquations, Graph, Digraph, Bigraph, Matrix)
(scoreEquations, Digraph, Bigraph, Graph, Matrix)
(scoreEquations, Bigraph, Graph, Digraph, Matrix)
(scoreEquations, Graph, Bigraph, Digraph, Matrix)
(scoreEquations, Bigraph, Digraph, Graph, Matrix)
(scoreEquations, Digraph, Graph, Bigraph, Matrix)
(scoreEquations, Matrix, MixedGraph)
(scoreEquations, Matrix, Graph)
(scoreEquations, Matrix, Digraph)
(scoreEquations, Matrix, Bigraph)
(scoreEquations, Matrix, Graph, Digraph)
(scoreEquations, Matrix, Digraph,Graph)
(scoreEquations, Matrix, Digraph,Bigraph)
(scoreEquations, Matrix, Bigraph,Digraph)
(scoreEquations, Matrix, Graph, Bigraph)
(scoreEquations, Matrix, Bigraph, Graph)
(scoreEquations, Matrix, Graph, Digraph, Bigraph)
(scoreEquations, Matrix, Digraph, Bigraph, Graph)
(scoreEquations, Matrix, Bigraph, Graph, Digraph)
(scoreEquations, Matrix, Graph, Bigraph, Digraph)
(scoreEquations, Matrix, Bigraph, Digraph, Graph)
(scoreEquations, Matrix, Digraph, Graph, Bigraph)
Headline
score equations of the log-likelihood function of a Gaussian graphical model
Usage
scoreEquations(R,U)
Inputs
R:Ring
defined as a @TO gaussianRing@ of @TO Graph@, or @TO Digraph@, or @TO Bigraph@, or @TO MixedGraph@
U:Matrix
or @TO List@ of sample data.
Alternatively, the input can be the sample covariance @TO Matrix@ by setting the optional input @TO SampleData@ to false
Outputs
:Sequence
consisting either of (Ideal) or (Ideal,Matrix)
where the ideal is generated by the score equations of the log-likelihood
function of the Gaussian model and the matrix is the symbolic covariance
matrix of the model
Description
Text
This function computes the score equations that arise from taking
partial derivatives of the log-likelihood function of the concentration matrix
(the inverse of the covariance matrix) of a Gaussian graphical
statistical model and returns the ideal generated by such equations.
The input of this function is a @TO gaussianRing@ and statistical data.
The latter can be given as a matrix or a list of observations. The rows of the matrix or the elements of the list are observation vectors given as lists.
It is possible to input the sample covariance matrix directly by using the optional input @TO SampleData@.
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
R = gaussianRing(G);
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
JU=scoreEquations(R,U)
V = sampleCovarianceMatrix U
JV=scoreEquations(R,V,SampleData=>false)
Text
@TO SaturateOptions@ allows to use all functionalities of @TO saturate@.
@TO Saturate@ determines whether to saturate. Note that the latter will not
provide the score equations of the model.
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
R = gaussianRing(G);
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
J=scoreEquations(R,U,SaturateOptions => {Strategy => Eliminate})
JnoSat=scoreEquations(R,U,Saturate=>false)
Text
The ML-degree of the model is the degree of the score equations ideal. The ML-degree
of the running example is 1:
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}});
R = gaussianRing(G);
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
J = scoreEquations(R,U)
dim J, degree J
///
doc ///
Key
CovarianceMatrix
Headline
optional input to output covariance matrix
SeeAlso
scoreEquations
///
doc ///
Key
[scoreEquations,CovarianceMatrix]
Headline
output covariance matrix
Usage
scoreEquations(R,U,CovarianceMatrix=>b)
Inputs
b:Boolean
default is false
Description
Text
@TO [scoreEquations,CovarianceMatrix]@ is set to false by default. If b is true,
@TO scoreEquations@ gives an additional output:
the covariance matrix with rational entries in the same variables as the ideal of score equations.
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
R=gaussianRing(G);
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
(J,Sigma)=scoreEquations(R,U,CovarianceMatrix=>true);
Sigma
SeeAlso
scoreEquations
///
doc ///
Key
Saturate
Headline
optional input whether to saturate
SeeAlso
scoreEquations
solverMLE
///
doc ///
Key
[scoreEquations, Saturate]
Headline
whether to saturate
Usage
scoreEquations(R,U,Saturate=>b)
Inputs
b:Boolean
default is true
Description
Text
@TO [scoreEquations,Saturate]@ is set to true by default. If b is false, saturation is not performed.
Avoiding saturation is only intended for big computations
when saturation cannot be computed or the computational time is very high.
When @TO Saturate@ is set to false, @TO scoreEquations@ might not output the ideal
generated by score equations because of the existence of vanishing denominators.
For graphs with only undirected edges, the ideal of score equations is the saturation of the
outputted ideal by the determinant of the concentration matrix. In the general case,
the ideal of score equations consist of the saturation of the outputted ideal by the product of denominators
of the Jacobian matrix.
For example, in the following case the degree of the ideal prior to saturation already gives the right ML-degree:
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
R=gaussianRing(G);
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
JnoSat=scoreEquations(R,U,Saturate=>false);
dim JnoSat
degree JnoSat
J=scoreEquations(R,U)
degree JnoSat==degree J
SeeAlso
scoreEquations
///
doc ///
Key
[solverMLE, Saturate]
Headline
whether to saturate
Usage
solverMLE(G,U,Saturate=>b)
Inputs
b:Boolean
default is true
Description
Text
@TO [solverMLE,Saturate]@ is set to true by default.
If we set @TO Saturate@ to false in @TO solverMLE@, saturation will not be
performed when computing the score equations of the log-likelihood function,
see @TO [scoreEquations,Saturate]@.
If the ideal returned by @TO scoreEquations@ has positive dimension, @TO solverMLE@
gives this ideal as output.
On the other hand, if we obtain a zero-dimensional ideal in @TO scoreEquations@,
@TO solverMLE@ computes the solutions of this polynomial system and returns:
* the maximal value of the log-likelihood function attained by positive definite matrices
corresponding to such solutions,
* the positive definite matrices where the maximum is attained,
* the degree of the zero-dimensional ideal.
Be aware that this output might not correspond to the right MLE.
Example
G=graph{{1,2},{2,3},{3,4},{1,4}}
U=random(ZZ^4,ZZ^4)
solverMLE(G,U,Saturate=>false)
SeeAlso
solverMLE
scoreEquations
[scoreEquations,Saturate]
///
doc ///
Key
SaturateOptions
Headline
optional input to use options "saturate"
SeeAlso
scoreEquations
solverMLE
Saturate
saturate
///
doc ///
Key
[scoreEquations, SaturateOptions]
Headline
use options from "saturate"
Usage
scoreEquations(R,U,SaturateOptions=>L)
Inputs
L: List
of options to set up saturation. Accepts any option from the function
@TO saturate@
Description
Text
Default @TO SaturateOptions@ in @TO scoreEquations@ are the
default options in @TO saturate@. All optional input in @TO saturate@ is allowed.
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}})
R=gaussianRing(G)
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
J=scoreEquations(R,U,SaturateOptions => {DegreeLimit=>1, MinimalGenerators => false})
SeeAlso
scoreEquations
Saturate
saturate
solverMLE
///
doc ///
Key
[solverMLE, SaturateOptions]
Headline
use options from "saturate"
Usage
solverMLE(G,U,SaturateOptions=>L)
Inputs
L: List
of options to set up saturation. Accepts any option from the function
@TO saturate@
Description
Text
Default @TO SaturateOptions@ in @TO solverMLE@ are the default options in @TO saturate@.
All optional input in @TO saturate@ is allowed.
Example
G=graph{{1,2},{2,3},{3,4},{1,4}}
U=random(ZZ^4,ZZ^4)
solverMLE(G,U,SaturateOptions => {DegreeLimit=>1, MinimalGenerators => false})
SeeAlso
scoreEquations
Saturate
saturate
solverMLE
///
doc ///
Key
Solver
Headline
optional input to choose numerical solver
SeeAlso
solverMLE
EigenSolver
NumericalAlgebraicGeometry
zeroDimSolve
solveSystem
///
doc ///
Key
[solverMLE, Solver]
Headline
choose numerical solver
Usage
solverMLE(G,U,Solver=>P)
Inputs
P: String
name of the corresponding package
Description
Text
This option allows to choose which numerical solver to use to estimate the critical
points. There are two options: "EigenSolver" or "NAG4M2" (Numerical Algebraic Geometry for Macaulay2).
The default and strongly recommended option is "EigenSolver", in which case
the function @TO zeroDimSolve@ is used.
If "NAG4M2" is chosen, then @TO solveSystem@ is used.
Example
G=mixedGraph(graph{{a,b}},digraph{{a,d}},bigraph{{c,d}})
U=matrix{{1, 2, 5, 1}, {5, 3, 2, 1}, {4, 3, 5, 10}, {2, 5,1, 3}};
solverMLE (G,U,Solver=>"EigenSolver")
solverMLE (G,U,Solver=>"NAG4M2")
SeeAlso
solverMLE
EigenSolver
NumericalAlgebraicGeometry
zeroDimSolve
solveSystem
///
doc ///
Key
OptionsEigenSolver
Headline
optional input to use options of "zeroDimSolve" in "EigenSolver"
SeeAlso
solverMLE
EigenSolver
zeroDimSolve
///
doc ///
Key
[solverMLE, OptionsEigenSolver]
Headline
use options of "zeroDimSolve" in "EigenSolver"
Usage
solverMLE(G,U,Solver=>"EigenSolver",OptionsEigenSolver=>L)
Inputs
L: List
of optional inputs in @TO zeroDimSolve@
Description
Text
Default @TO OptionsEigenSolver@ in @TO solverMLE@ are the default options in @TO zeroDimSolve@ in @TO EigenSolver@.
All optional input in @TO zeroDimSolve@ is allowed.
Example
G=mixedGraph(graph{{a,b},{b,c}})
solverMLE(G,matrix{{1,0,0},{0,1,0},{0,0,1}},Solver=>"EigenSolver",OptionsEigenSolver=>{Multiplier =>1, Strategy=>"Stickelberger"})
Text
In fact, since @TO zeroDimSolve@ is the current default numerical solver.
the sintaxis below is enough:
Example
G=mixedGraph(graph{{a,b},{b,c}})
solverMLE(G,matrix{{1,0,0},{0,1,0},{0,0,1}},OptionsEigenSolver=>{Multiplier =>1, Strategy=>"Stickelberger"})
SeeAlso
solverMLE
EigenSolver
zeroDimSolve
///
doc ///
Key
OptionsNAG4M2
Headline
optional parameter to set the parameters of solveSystem in NumericalAlgebraicGeometry
SeeAlso
solverMLE
NumericalAlgebraicGeometry
solveSystem
///
doc ///
Key
[solverMLE, OptionsNAG4M2]
Headline
use options of "solveSystem" in "NumericalAlgebraicGeometry"
Usage
solverMLE(G,U,Solver=>"NAG4M2",OptionsNAG4M2=>L)
Inputs
L: List
of optional inputs to @TO solveSystem@
Description
Text
Default @TO OptionsNAG4M2@ in @TO solverMLE@ when setting @TO [solverMLE,Solver]@
to "NAG4M2" are the default options of @TO solveSystem@ in @TO NumericalAlgebraicGeometry@.
All optional input in @TO solveSystem@ is allowed.
Example
G=mixedGraph(graph{{a,b},{b,c}})
solverMLE(G,matrix{{1,0,0},{0,1,0},{0,0,1}},Solver=>"NAG4M2",OptionsNAG4M2=>{tStep =>.01,numberSuccessesBeforeIncrease => 5})
SeeAlso
solverMLE
NumericalAlgebraicGeometry
solveSystem
///
doc ///
Key
RealPrecision
Headline
optional input to choose the number of decimals used to round to QQ when inputing data in RR
SeeAlso
scoreEquations
solverMLE
///
doc ///
Key
[scoreEquations, RealPrecision]
Headline
the number of bits of precision to use in the computation
Usage
scoreEquations(R,U,RealPrecision=>n)
Inputs
n: ZZ
default is 53
Description
Text
This optional input only applies when the sample data or the sample covariance matrix has real entries.
By default, the precision is 53 (the default precision for real numbers in M2).
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
R=gaussianRing(G);
U = matrix{{6.2849049, 10.292875, 1.038475, 1.1845757}, {3.1938475, 3.2573, 1.13847, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
J=scoreEquations(R,U,RealPrecision=>3)
SeeAlso
scoreEquations
///
doc ///
Key
[solverMLE, RealPrecision]
Headline
the number of bits of precision to use in the computation
Usage
solverMLE(G,U,RealPrecision=>n)
Inputs
n: ZZ
default is 53
Description
Text
This optional input only applies when the sample data or the sample covariance matrix has real entries.
By default, the precision is 53 (the default precision for real numbers in M2).
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
U = matrix{{6.2849049, 10.292875, 1.038475, 1.1845757}, {3.1938475, 3.2573, 1.13847, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}};
solverMLE(G,U,RealPrecision=>10)
SeeAlso
scoreEquations
solverMLE
///
doc ///
Key
SampleData
Headline
optional input to allow to input the sample covariance matrix instead of sample data
SeeAlso
scoreEquations
solverMLE
///
doc ///
Key
[scoreEquations, SampleData]
Headline
input sample covariance matrix instead of sample data
Usage
scoreEquations(R,U,SampleData=>b)
Inputs
b: Boolean
default is true
Description
Text
@TO scoreEquations@ requires a matrix or a list of sample data as part of the default input.
Setting @TO [scoreEquations,SampleData]@ to false allows the user to enter a sample
covariance matrix instead of sample data. It must be a symmetric matrix.
Example
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
R=gaussianRing(G);
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}}
J=scoreEquations(R,U,SampleData=>true)
V=sampleCovarianceMatrix(U)
I=scoreEquations(R,V,SampleData=>false)
SeeAlso
scoreEquations
///
doc ///
Key
[solverMLE, SampleData]
Headline
input sample covariance matrix instead of sample data
Usage
solverMLE(G,U,SampleData=>b)
Inputs
b: Boolean
default is true
Description
Text
@TO solverMLE@ requires a matrix or a list of sample data as part of the default input.
Setting @TO [solverMLE,SampleData]@ to false allows the user to enter a sample
covariance matrix as input. It must be a symmetric matrix.
Example
G=graph{{1,2},{2,3},{3,4},{1,4}}
U=random(ZZ^4,ZZ^4)
solverMLE(G,U,SampleData=>true)
V=sampleCovarianceMatrix(U)
solverMLE(G,V,SampleData=>false)
SeeAlso
scoreEquations
solverMLE
///
doc ///
Key
ConcentrationMatrix
Headline
optional input to output MLE for concentration matrix instead of MLE for covariance matrix
SeeAlso
solverMLE
///
doc ///
Key
[solverMLE, ConcentrationMatrix]
Headline
output MLE for concentration matrix instead of MLE for covariance matrix
Usage
solverMLE(G,U,ConcentrationMatrix=>b)
Inputs
b: Boolean
false by default
Description
Text
By default, @TO solverMLE@ outputs the MLE for the covariance matrix. By providing true as the value for the option
@TO [solverMLE, ConcentrationMatrix]@, @TO solverMLE@ provides the MLE for
the concentration matrix.
Note that both the maximum value attained in the log-likelihood function and
the ML-degree remain the same.
Example
G= mixedGraph(graph{{a,b},{b,c}},digraph {{a,d},{c,e},{f,g}},bigraph {{d,e}})
solverMLE (G, random(QQ^7,QQ^7))
solverMLE (G, random(QQ^7,QQ^7), ConcentrationMatrix=>true)
SeeAlso
solverMLE
///
doc ///
Key
checkPD
(checkPD,List)
(checkPD,Matrix)
Headline
returns positive definite matrices from a list of symmetric matrices
Usage
checkPD(L)
Inputs
L: List
of symmetric matrices, or a single symmetric @TO Matrix@
Outputs
: List
of positive definite matrices
Description
Text
This function takes a list of symmetric matrices (or a single symmetric matrix) and returns the sublist of its positive definite matrices.
If there are no positive definite matrices in the list, it returns an empty list.
Note that the function does not check whether the matrices in the original list are symmetric.
If a matrix contains an imaginary part below the tolerance level, then only
the real part is reported in the output. (See @TO [checkPD, ZeroTolerance]@)
Example
L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}}}
checkPD(L)
///
doc ///
Key
checkPSD
(checkPSD,List)
(checkPSD,Matrix)
Headline
returns positive semidefinite matrices from a list of symmetric matrices
Usage
checkPSD(L)
Inputs
L: List
of symmetric matrices, or a single symmetric @TO Matrix@
Outputs
: List
of positive semidefinite matrices
Description
Text
This function takes a list of symmetric matrices (or a single symmetric matrix) and returns the sublist of its
positive semidefinite matrices.
Note that the function does not check whether the matrices in the original list are symmetric.
If there are no positive semidefinite matrices in the list, it returns an empty list.
Example
L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}},matrix{{0,0},{0,0}}}
checkPSD(L)
///
doc ///
Key
ZeroTolerance
Headline
optional input to set the largest absolute value that should be treated as zero
SeeAlso
checkPD
checkPSD
solverMLE
///
doc ///
Key
[solverMLE, ZeroTolerance]
Headline
optional input to set the largest absolute value that should be treated as zero
Usage
solverMLE(G,U,ZeroTolerance=>n)
Inputs
n: RR
default is 1e-10
Description
Text
After computing the variety of the zero-dimensional ideal of the score equations,
{\tt solverMLE} needs to determine which points lie in the PD cone. A matrix is
assumed to be positive definite if for all eigenvalues e:
- @TO realPart@ e > ZeroTolerance
- @TO abs @ @TO imaginaryPart@ e <= ZeroTolerance
If an MLE matrix contains an imaginary part below the tolerance level, then only
the real part is reported in the output. (See @TO [checkPD, ZeroTolerance]@)
SeeAlso
checkPD
///
doc ///
Key
[checkPD, ZeroTolerance]
Headline
optional input to set the largest absolute value that should be treated as zero
Usage
checkPD(L,ZeroTolerance=>n)
Inputs
n: RR
default is 1e-10
Description
Text
A matrix is assumed to be positive definite if for all eigenvalues e:
- @TO realPart@ e > ZeroTolerance
- @TO abs @ @TO imaginaryPart@ e <= ZeroTolerance
If a matrix contains an imaginary part below the tolerance level, then only
the real part is reported in the output.
Example
L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}},
matrix{{10^(-10)+10^(-10)*sqrt(-1),0},{0,10^(-10)+10^(-10)*sqrt (-1)}},
matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},
matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}
}
checkPD L
SeeAlso
[checkPSD, ZeroTolerance]
solverMLE
///
doc ///
Key
[checkPSD, ZeroTolerance]
Headline
optional input to set the largest absolute value that should be treated as zero
Usage
checkPD(L,ZeroTolerance=>n)
Inputs
n: RR
default is 1e-10
Description
Text
A matrix is assumed to be positive semidefinite if for all eigenvalues e:
- @TO realPart@ e >= -ZeroTolerance
- @TO abs @ @TO imaginaryPart@ e <= ZeroTolerance
If a matrix contains an imaginary part below the tolerance level, then only
the real part is reported in the output.
Example
L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}},
matrix{{10^(-10)+10^(-10)*sqrt(-1),0},{0,10^(-10)+10^(-10)*sqrt (-1)}},
matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},
matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}
}
checkPD L
SeeAlso
[checkPD, ZeroTolerance]
///
doc ///
Key
MLdegree
(MLdegree,Ring)
Headline
ML-degree of a graphical model
Usage
MLdegree(R)
Inputs
R: Ring
defined as a @TO gaussianRing@ of @TO Graph@, or @TO Digraph@, or @TO Bigraph@, or @TO MixedGraph@
Outputs
: ZZ
the ML-degree of the model
Description
Text
This function computes the ML-degree of a graphical model. It takes as input
a @TO gaussianRing@ of a @TO Graph@, or a @TO Digraph@, or a @TO Bigraph@, or a @TO MixedGraph@.
It computes the degree of the score equation ideal given by @TO scoreEquations@
with a random sample data matrix.
We compute the ML-degree of the 4-cycle:
Example
G=graph{{1,2},{2,3},{3,4},{4,1}}
MLdegree(gaussianRing G)
///
doc ///
Key
solverMLE
(solverMLE, Ring, List)
(solverMLE, Ring, Matrix)
(solverMLE, List, Ring)
(solverMLE, Matrix, Ring)
(solverMLE,MixedGraph,List)
(solverMLE, MixedGraph, Matrix)
(solverMLE, Graph, List)
(solverMLE, Digraph, List)
(solverMLE, Bigraph, List)
(solverMLE, Graph, Digraph,List)
(solverMLE, Digraph, Graph,List)
(solverMLE, Digraph, Bigraph, List)
(solverMLE, Bigraph, Digraph, List)
(solverMLE, Graph, Bigraph, List)
(solverMLE, Bigraph, Graph, List)
(solverMLE, Graph, Digraph, Bigraph, List)
(solverMLE, Digraph, Bigraph, Graph, List)
(solverMLE, Bigraph, Graph, Digraph, List)
(solverMLE, Graph, Bigraph, Digraph, List)
(solverMLE, Bigraph, Digraph, Graph, List)
(solverMLE, Digraph, Graph, Bigraph, List)
(solverMLE, List, MixedGraph)
(solverMLE, List, Graph)
(solverMLE, List, Digraph)
(solverMLE, List, Bigraph)
(solverMLE, List, Graph, Digraph)
(solverMLE, List, Digraph,Graph)
(solverMLE, List, Digraph,Bigraph)
(solverMLE, List, Bigraph,Digraph)
(solverMLE, List, Graph, Bigraph)
(solverMLE, List, Bigraph, Graph)
(solverMLE, List, Graph, Digraph, Bigraph)
(solverMLE, List, Digraph, Bigraph, Graph)
(solverMLE, List, Bigraph, Graph, Digraph)
(solverMLE, List, Graph, Bigraph, Digraph)
(solverMLE, List, Bigraph, Digraph, Graph)
(solverMLE, List, Digraph, Graph, Bigraph)
(solverMLE, Graph, Matrix)
(solverMLE, Digraph, Matrix)
(solverMLE, Bigraph, Matrix)
(solverMLE, Graph, Digraph,Matrix)
(solverMLE, Digraph, Graph,Matrix)
(solverMLE, Digraph, Bigraph, Matrix)
(solverMLE, Bigraph, Digraph, Matrix)
(solverMLE, Graph, Bigraph, Matrix)
(solverMLE, Bigraph, Graph, Matrix)
(solverMLE, Graph, Digraph, Bigraph, Matrix)
(solverMLE, Digraph, Bigraph, Graph, Matrix)
(solverMLE, Bigraph, Graph, Digraph, Matrix)
(solverMLE, Graph, Bigraph, Digraph, Matrix)
(solverMLE, Bigraph, Digraph, Graph, Matrix)
(solverMLE, Digraph, Graph, Bigraph, Matrix)
(solverMLE, Matrix, MixedGraph)
(solverMLE, Matrix, Graph)
(solverMLE, Matrix, Digraph)
(solverMLE, Matrix, Bigraph)
(solverMLE, Matrix, Graph, Digraph)
(solverMLE, Matrix, Digraph,Graph)
(solverMLE, Matrix, Digraph,Bigraph)
(solverMLE, Matrix, Bigraph,Digraph)
(solverMLE, Matrix, Graph, Bigraph)
(solverMLE, Matrix, Bigraph, Graph)
(solverMLE, Matrix, Graph, Digraph, Bigraph)
(solverMLE, Matrix, Digraph, Bigraph, Graph)
(solverMLE, Matrix, Bigraph, Graph, Digraph)
(solverMLE, Matrix, Graph, Bigraph, Digraph)
(solverMLE, Matrix, Bigraph, Digraph, Graph)
(solverMLE, Matrix, Digraph, Graph, Bigraph)
Headline
Maximum likelihood estimate of a loopless mixed graph
Usage
solverMLE(G,U)
Inputs
G:Graph
or @ofClass Digraph@, or @ofClass Bigraph@, or @ofClass MixedGraph@
U:Matrix
or @ofClass List@ of sample data.
Alternatively, the sample covariance matrix can be given as input
by setting @TO SampleData@ to false
Outputs
: Sequence
of length 3, whose first element is the maximum value attained in the log-likelihood function (of type @TO RR@),
its second element is the MLE (or MLEs) of the covariance matrix (of types @TO Matrix@ or @TO List@),
and its third element is the ML-degree of the model (of type @TO ZZ@).
By providing true as the value for the option @TO ConcentrationMatrix@, the MLE for the concentration matrix can be given as output.
Description
Text
This function takes as input a @TO2{Graph,"graph"}@, or a @TO2{Digraph,"digraph"}@, or a @TO2{Bigraph,"bigraph"}@ or a @TO2{MixedGraph,"mixed graph"}@ and a list or matrix that encodes, by default, the sample data.
It computes the critical points of the score equations and
selects the maximum value achieved among those that lie in the cone of positive-definite matrices.
The default output is the maximum value in the log-likelihood function, maximum likelihood estimators (MLE) for the covariance matrix
and the ML-degree of the model.
MLE for the concentration matrix can be obtained by setting the optional input @TO ConcentrationMatrix@ to false.
The same optional inputs as in @TO scoreEquations@ can be used, plus extra optional inputs related to
the numerical solver (EigenSolver by default, NAG4M2 alternatively) and its functionalities.
Below we reproduce Example 2.1.13 for the 4-cycle in the book: Mathias Drton, Bernd Sturmfels and Seth Sullivant,
{\em Lectures on Algebraic Statistics}, Oberwolfach Seminars, Vol 40, Birkhauser, Basel, 2009.
Example
G=graph{{1,2},{2,3},{3,4},{1,4}};
U =matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
solverMLE(G,U)
Text
The data sample can also be given as a list:
Example
G=graph{{1,2},{2,3},{3,4},{1,4}};
U = {{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
solverMLE(G,U)
Text
In the following example we compute the MLE for the covariance matrix of
the graphical model associated to the graph $1\rightarrow 2,1\rightarrow 3,2\rightarrow 3,3\rightarrow 4,3<-> 4$
In this case we give as input the sample covariance matrix:
Example
G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}});
S = matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}};
solverMLE(G,S,SampleData=>false)
Text
Next we provide the MLE for the concentration matrix of the graphical model
associated to the graph $1\rightarrow 3,2\rightarrow 4,3<->4,1 - 2$.
Again the sample covariance matrix is given as input.
Example
G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}},graph {{1,2}});
S = matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}};
solverMLE(G,S,SampleData=>false,ConcentrationMatrix=>true)
Text
{\bf Application to positive definite matrix completion problems}
Text
Consider the following symmetric matrix with some unknown entries:
Example
R=QQ[x,y];
M=matrix{{115,-13,x,47},{-13,5,7,y},{x,7,27,-21},{47,y,-21,29}}
Text
Unknown entries correspond to non-edges of the 4-cycle. A positive definite completion of this matrix
is obtained by giving values to x and y and computing the MLE for the covariance matrix in the Gaussian graphical model
given by the 4-cycle. To understand which values of x and y will result in a maximum likelihood estimate,
see Example 12.16 in the book: Mateusz Michalek and Bernd Sturmfels,
{\em Invitation to Nonlinear Algebra}, Graduate Studies in Mathematics, Vol ???, American Mathematical Society, 2021.
Example
G=graph{{1,2},{2,3},{3,4},{1,4}};
V=matrix{{115,-13,-29,47},{-13,5,7,-11},{-29,7,27,-21},{47,-11,-21,29}}
(mx,MLE,ML)=solverMLE(G,V,SampleData=>false)
Text
The MLE of the covariance matrix is the unique positive definite completion of the matrix M such that its inverse, namely the concentration matrix,
has zero's in the entries corresponding to non-edges of the graph.
Observe that all entries of V remain the same in the MLE except for those that correspond to non-edges of the graph.
SeeAlso
checkPD
checkPSD
scoreEquations
jacobianMatrixOfRationalFunction
sampleCovarianceMatrix
///
--******************************************--
-- TESTS --
--******************************************--
TEST /// --test jacobianMatrixOfRationalFunction
R=QQ[x,y];
FR=frac R;
F=1/(x^2+y^2);
M=entries jacobianMatrixOfRationalFunction(F);
N=transpose {{-2*x/(x^2 + y^2)^2,-2*y/(x^2 + y^2)^2 }};
assert(M === N)
///
TEST /// --test jacobianMatrixOfRationalFunction
R=QQ[x_1,x_2,x_3];
FR=frac R;
M=entries jacobianMatrixOfRationalFunction( (x_1^2*x_2)/(x_1+x_2^2+x_3^3) );
N=transpose {{2*x_1*x_2/(x_2^2 + x_3^3 + x_1) - x_1^2*x_2/(x_2^2 + x_3^3 + x_1)^2, -2*x_1^2*x_2^2/(x_2^2 + x_3^3 + x_1)^2 + x_1^2/(x_2^2 + x_3^3 + x_1) , -3*x_1^2*x_2*x_3^2/(x_2^2 + x_3^3 + x_1)^2 }};
assert(M === N)
///
TEST /// --test sampleCovarianceMatrix
M = matrix{{1, 2, 0},{-1, 0, 5/1},{3, 5, 2/1},{-1, -4, 1/1}};
N = sampleCovarianceMatrix(M);
A = matrix {{11/4, 39/8, -1}, {39/8, 171/16, 0}, {-1, 0, 7/2}};
assert(N===A)
///
TEST /// --test sampleCovarianceMatrix with sample of bigger size than the covariance matrix
X = matrix{{36, -3, -25, -36},{-10, 11, -29, -20},{-18, 33, -15, -11},{-42, 0, 20, 43}, {-30, -26, 32, 2},{2, -38, -24, -43}};
Y = sampleCovarianceMatrix(X);
B = matrix {{5621/9, -1037/18, -7835/18, -10565/18}, {-1037/18, 19505/36, -4897/36, 5147/36}, {-7835/18, -4897/36, 20465/36, 18941/36}, {-10565/18, 5147/36, 18941/36, 28889/36}};
assert(Y===B)
///
TEST /// --test sampleCovarianceMatrix with list input
X = {{48,89,27,28},{23,19,29,94},{135,23,44,71},{91,75,24,98}};
Y = sampleCovarianceMatrix(X);
B = matrix {{29147/16, -1313/8, 220, 1609/16}, {-1313/8, 3827/4, -155, -3451/8}, {220, -155, 119/2, -63/4}, {1609/16, -3451/8, -63/4, 12379/16}};
assert(Y===B)
///
TEST /// --test scoreEquations for a mixedGraph with directed and bidirected edges
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph {{3,4}});
R=gaussianRing(G);
U = matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}};
J=scoreEquations(R,U);
I=ideal(20*p_(3,4)+39,50*p_(4,4)-271,440104*p_(3,3)-742363,230*p_(2,2)-203,16*p_(1,1)-115,5*l_(3,4)+2,110026*l_(2,3)-2575,55013*l_(1,3)-600,115*l_(1,2)+26);
assert(J===I)
///
TEST /// --test scoreEquations for a graph (and random input data)
G=graph{{1,2},{2,3},{3,4},{1,4}};
R=gaussianRing(G);
U=random(ZZ^4,ZZ^4);
J=scoreEquations(R,U);
assert(dim J===0)
assert(degree J===5)
///
TEST /// --score equations of sample data equals score equations of its sample covariance data with SampleData=>false
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}})
R = gaussianRing(G)
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}}
JU=scoreEquations(R,U)
RU=ring(JU)
V = sampleCovarianceMatrix U
JV=scoreEquations(R,V,SampleData=>false)
JV=sub(JV,RU)
assert(JU==JV)
///
TEST /// --score equations with elimination strategy equals default saturation strategy
G = mixedGraph(digraph {{1,2},{1,3},{2,3},{3,4}},bigraph{{3,4}})
R = gaussianRing(G)
U = matrix{{6, 10, 1/3, 1}, {3/5, 3, 1/2, 1}, {4/5, 3/2, 9/8, 3/10}, {10/7, 2/3,1, 8/3}}
JU=scoreEquations(R,U)
RU=ring(JU)
J=scoreEquations(R,U,SaturateOptions => {Strategy => Eliminate})
J=sub(J,RU)
assert(J==JU)
///
TEST /// --test checkPD
L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}},matrix{{0.0001*sqrt(-1),0},{0,0.0000001*sqrt (-1)}}};
Y = checkPD(L);
B = {matrix{{1, 0}, {0, 1}}};
assert(Y===B)
///
TEST /// --test ZeroTolerance checkPD
L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}}, matrix{{10^(-10)+10^(-10)*sqrt(-1),0},{0,10^(-10)+10^(-10)*sqrt (-1)}},matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}}
assert(checkPD L ==={matrix {{1e-9, 0}, {0, 1e-9}}, sub(matrix {{1, 0}, {0, 1}},RR)})
///
TEST /// --test ZeroTolerance checkPSD
L={matrix{{10^(-9)+10^(-10)*sqrt(-1),0},{0,10^(-9)+10^(-10)*sqrt (-1)}}, matrix{{-10^(-10)+10^(-10)*sqrt(-1),0},{0,-10^(-10)+10^(-10)*sqrt (-1)}},matrix{{1+10^(-10)*sqrt(-1),0},{0,1+10^(-10)*sqrt (-1)}},matrix{{1-10^(-9)*sqrt(-1),0},{0,1+10^(-9)*sqrt (-1)}}}
assert(checkPSD L ==={matrix {{1e-9, 0}, {0, 1e-9}},matrix {{-1e-10, 0}, {0, -1e-10}}, sub(matrix {{1, 0}, {0, 1}},RR)})
///
TEST /// --test checkPSD
L={matrix{{1,0},{0,1}},matrix{{-2,0},{0,1}},matrix{{sqrt(-1),0},{0,sqrt (-1)}},matrix{{0.0001*sqrt(-1),0},{0,0.0000001*sqrt (-1)}},matrix{{0,0},{0,0}}};
Y = checkPSD(L);
B = {matrix{{1, 0}, {0, 1}},matrix{{0,0},{0,0}}};
assert(Y===B)
///
TEST /// --test MLdegree
G=graph{{1,2},{2,3},{3,4},{4,1}}
assert(MLdegree(gaussianRing G)==5)
///
TEST /// --test solverMLE graph
G=graph{{1,2},{2,3},{3,4},{1,4}}
U =matrix{{1,2,1,-1},{2,1,3,0},{-1, 0, 1, 1},{-5, 3, 4, -6}}
(mx,MLE,ML)=solverMLE(G,U)
assert(round(4,mx)==-6.2615)
assert(ML==5)
///
TEST /// --test solverMLE for mixedGraph with directed and bidirected edges
G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}})
S = matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}}
(mx,MLE,ML)=solverMLE(G,S,SampleData=>false)
assert(ML==5)
assert(round(5, mx)==-1.14351)
///
TEST /// --test solverMLE for mixedGraph with all kind of edges and concentration matrix
G = mixedGraph(digraph {{1,3},{2,4}},bigraph {{3,4}},graph {{1,2}})
S = matrix {{7/20, 13/50, -3/50, -19/100}, {13/50, 73/100, -7/100, -9/100},{-3/50, -7/100, 2/5, 3/50}, {-19/100, -9/100, 3/50, 59/100}}
(mx,MLE,ML)=solverMLE(G,S,SampleData=>false,ConcentrationMatrix=>true)
assert(ML==5)
assert(round(4,mx)==-.8362)
///
TEST /// --test solverMLE graph with complete test
G=graph{{1,2},{2,3},{3,4},{1,4}}
V =matrix{{5,1,3,2},{1,5,1,6},{3,1,5,1},{2,6,1,5}}
(mx,MLE,ML)=solverMLE(G,V,SampleData=>false,ConcentrationMatrix=>true)
assert(round(4,mx)==-10.1467)
assert(ML==5)
assert(round(6,MLE_(0,1))== -.038900)
assert(round(6,MLE_(3,1))== 0)
///
TEST /// -- test solverMLE with NAG4M2
G=mixedGraph(graph{{a,b}},digraph{{a,d}},bigraph{{c,d}})
U=matrix{{1, 2, 5, 1}, {5, 3, 2, 1}, {4, 3, 5, 10}, {2, 5,1, 3}}
(mx,MLE,ML)= solverMLE (G,U,Solver=>"NAG4M2")
assert(round(5,mx)==-8.4691)
assert(ML==1)
assert(round(6,MLE_(2,0))==0)
assert(round(6,MLE_(1,1))== 1.1875)
assert(round(6,MLE_(3,2))== 3.264929)
///
--------------------------------------
--------------------------------------
end--
--------------------------------------
--------------------------------------