=========
Maximum Likelihood Degree
=========

Below is the computation (in Macaulay2) for ML-degree for linear concentration models. After reading the code, the function phi(a,n) computes the ML-degree for a general a dimensional subspace of n x n symmetric matrices.


loadPackage "SchurRings"
num=(n,j,k)->(
di=(n+1)*n/2-1;
pie=toList j_0;
ost=toList j_(#j-1);
if (#pie)==1 then H1=2*2^(pie_0) else H1=2^(pie_0)+2^(pie_1);
if (#ost)==1 then H2=-2*2^(ost_0) else H2=-2^(ost_0)-2^(ost_1);
return (-1)^di*H1^k*H2^(di-k);
);

sas=(s1,s2)->(
s1=toList s1;
s2=toList s2;
if ((#s1==1) and (#s2==1)) then return (-2*2^(s1_0)+2*2^(s2_0));
if ((#s1==1) and (#s2==2)) then return (-2*2^(s1_0)+2^(s2_1)+2^(s2_0));
if ((#s1==2) and (#s2==1)) then return (2*2^(s2_0)-2^(s1_1)-2^(s1_0));
if ((#s1==2) and (#s2==2)) then return (-2^(s1_0)-2^(s1_1)+2^(s2_0)+2^(s2_1));
);

dwu=j->(
d=1;
dl=#j;
for i from 0 to (dl-1) do(
jj=toList j_i;
if (#(jj)==2) then d=d*(2^(jj_0)-2^(jj_1))*(-2^(jj_0)+2^(jj_1));
);
return d;
);

denom=(n,j)->(
dl=#j;
d=1;
for i from 0 to (dl-2) do(
 d=d*sas(j_i,j_(i+1));
);
ntka=set toList(1..n);
for i in j do (
 ntka=ntka-i;
 for pier in (toList i) do ( 
  for dr in (toList ntka) do (
   d=d*(-2^pier+2^dr)
  );
 );
);
return d*dwu(j);
);

phi=(k,n)->(
 q=0;
 PAR=partitions(n,2);
 for p in PAR do (
  p=toList p;
  zzst=partitions(set toList (1..n),p);
  for i in zzst do (
   zz=permutations(toList i);
   for j in zz do (
    q=q+num(n,j,k)/denom(n,j);
   );
  );
 );
return q;
);