--Gaussian likelihood
restart
R = QQ[k11,k12,k22,k23,k24,k33,k34,k44]
K = matrix {{k11,k12,0,0},{k12,k22,k23,k24},{0,k23,k33,k34},{0,k24,k34,k44}}
X = matrix for i to 3 list for j to 3 list random(30)
S = X*transpose(X)
M1 = jacobian(ideal(det(K)));
M2 = det(K)*jacobian(ideal(trace(K*S)));
I = ideal (M1-M2);
J = saturate(I,det(K))


--Gaussian likelihood 2
restart
R = QQ[k11,k12,k22,k23,k24,k33,k34,k44,s11,s12,s13,s14,s22,s23,s24,s33,s34,s44]
K = matrix {{k11,k12,0,0},{k12,k22,k23,k24},{0,k23,k33,k34},{0,k24,k34,k44}}
Sigma = matrix {{s11,s12,s13,s14},{s12,s22,s23,s24},{s13,s23,s33,s34},{s14,s24,s34,s44}}
I1 = ideal (K*Sigma - identity(1))
X = matrix for i to 3 list for j to 3 list random(30)
S = X*transpose(X)
I2 = ideal(Sigma_(0,0)-S_(0,0),Sigma_(0,1)-S_(0,1),Sigma_(1,1)-S_(1,1),Sigma_(1,2)-S_(1,2),Sigma_(1,3)-S_(1,3),Sigma_(2,2)-S_(2,2),Sigma_(2,3)-S_(2,3),Sigma_(3,3)-S_(3,3))
I = I1 + I2
J = eliminate(I,{k11,k12,k22,k23,k24,k33,k34,k44})

--discrete MLE
restart
R1 = QQ[p_(0,0,0,0)..p_(1,1,1,1)]
R2 = QQ[p_(0,0,0,0)..p_(1,1,1,1),a_(0,0)..a_(1,1),b_(0,0,0)..b_(1,1,1)]
IF = ideal flatten flatten flatten for i to 1 list for j to 1 list for k to 1 list for l to 1 list p_(i,j,k,l)-a_(i,j)*b_(j,k,l)
JF = eliminate(IF,join(toList(a_(0,0)..a_(1,1)), toList(b_(0,0,0)..b_(1,1,1))))
JF = sub(JF,R1)
use R1
A = matrix{{ 1,1,1,1,0,0,0,0,0,0,0,0,0,0,0,0},
{0,0,0,0,1,1,1,1,0,0,0,0,0,0,0,0},
{0,0,0,0,0,0,0,0,1,1,1,1,0,0,0,0},
{0,0,0,0,0,0,0,0,0,0,0,0,1,1,1,1},
{1,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0},
{0,1,0,0,0,0,0,0,0,1,0,0,0,0,0,0},
{0,0,1,0,0,0,0,0,0,0,1,0,0,0,0,0},
{0,0,0,1,0,0,0,0,0,0,0,1,0,0,0,0},
{0,0,0,0,1,0,0,0,0,0,0,0,1,0,0,0},
{0,0,0,0,0,1,0,0,0,0,0,0,0,1,0,0},
{0,0,0,0,0,0,1,0,0,0,0,0,0,0,1,0},
{0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,1}}
U = transpose matrix {for i to 15 list random(30)}
P = transpose matrix {toList(p_(0,0,0,0)..p_(1,1,1,1))}
I = JF + ideal (A*U-A*P)
dim I
degree I