>
 

MATH 392 -- Seminar in Computational Commutative Algebra 

 

Maximum Likelihood Estimation in a Linear Mixture Model 

 

We wish to determine the values of theta[1] and theta[2]that maximize 

the log-likelihood function  

 

  l = log(p[A]^u[A]*p[C]^u[C]*p[G]^u[G]*p[T]^u[T]) 

 

where 

 

p[A] = .15*theta[1]+.27*theta[2]+`*`(.25, 1-theta[1]-theta[2]) 

p[C] = .33*theta[1]+.24*theta[2]+`*`(.25, 1-theta[1]-theta[2]) 

p[G] = .36*theta[1]+.23*theta[2]+`*`(.25, 1-theta[1]-theta[2])*p[T] and .36*theta[1]+.23*theta[2]+`*`(.25, 1-theta[1]-theta[2])*p[T] = .16*theta[1]+.26*theta[2]+`*`(.25, 1-theta[1]-theta[2])
p[G] = .36*theta[1]+.23*theta[2]+`*`(.25, 1-theta[1]-theta[2])*p[T] and .36*theta[1]+.23*theta[2]+`*`(.25, 1-theta[1]-theta[2])*p[T] = .16*theta[1]+.26*theta[2]+`*`(.25, 1-theta[1]-theta[2]) 

 

 

and from the data  u[A] = 10, u[C] = 14, u[G] = 15, u[T] = 10. 

> p[A] := `*`(15, 1/100)*theta[1]+27/100*theta[2]+`*`(25, 1/100)*(1-theta[1]-theta[2]); 1
 

-1/10*theta[1]+1/50*theta[2]+1/4 

> p[C] := 33/100*theta[1]+`*`(24, 1/100)*theta[2]+`*`(25, 1/100)*(1-theta[1]-theta[2]); 1
 

2/25*theta[1]-1/100*theta[2]+1/4 

> p[G] := `*`(36, 1/100)*theta[1]+23/100*theta[2]+`*`(25, 1/100)*(1-theta[1]-theta[2]); 1
 

11/100*theta[1]-1/50*theta[2]+1/4 

> p[T] := `*`(16, 1/100)*theta[1]+`*`(26, 1/100)*theta[2]+`*`(25, 1/100)*(1-theta[1]-theta[2]); 1
 

-9/100*theta[1]+1/100*theta[2]+1/4 

The log-likelihood function is 

> l := ln(p[A]^10*p[C]^14*p[G]^15*p[T]^10); 1
 

ln((-1/10*theta[1]+1/50*theta[2]+1/4)^10*(2/25*theta[1]-1/100*theta[2]+1/4)^14*(11/100*theta[1]-1/50*theta[2]+1/4)^15*(-9/100*theta[1]+1/100*theta[2]+1/4)^10)
ln((-1/10*theta[1]+1/50*theta[2]+1/4)^10*(2/25*theta[1]-1/100*theta[2]+1/4)^14*(11/100*theta[1]-1/50*theta[2]+1/4)^15*(-9/100*theta[1]+1/100*theta[2]+1/4)^10) 

To maximize, we apply the usual process from multivariable calculus: 

> eq1 := numer(simplify(diff(l, theta[1]))); 1
 

388080*theta[1]^3-1338*theta[2]^3-179706*theta[1]^2*theta[2]-2825625*theta[1]+457500*theta[2]-151950*theta[1]^2+42525*theta[1]*theta[2]-3825*theta[2]^2+27156*theta[1]*theta[2]^2+1359375
388080*theta[1]^3-1338*theta[2]^3-179706*theta[1]^2*theta[2]-2825625*theta[1]+457500*theta[2]-151950*theta[1]^2+42525*theta[1]*theta[2]-3825*theta[2]^2+27156*theta[1]*theta[2]^2+1359375 

> eq2 := numer(simplify(diff(l, theta[2]))); 1
 

-60100*theta[1]^3+196*theta[2]^3+27332*theta[1]^2*theta[2]+457500*theta[1]-77500*theta[2]+21150*theta[1]^2-6550*theta[1]*theta[2]+650*theta[2]^2-4052*theta[1]*theta[2]^2-218750
-60100*theta[1]^3+196*theta[2]^3+27332*theta[1]^2*theta[2]+457500*theta[1]-77500*theta[2]+21150*theta[1]^2-6550*theta[1]*theta[2]+650*theta[2]^2-4052*theta[1]*theta[2]^2-218750 

To solve the equationswe will use our Grobner basis tools: 

> with(Groebner); -1
 

> B := Basis([eq1, eq2], plex(theta[1], theta[2])); 1
 

[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]...
[15059072*theta[2]^9-204388812773437500*theta[2]^4-302623360342382812500*theta[2]^3-3313234506368408203125*theta[2]^2+14248392817565917968750*theta[2]-2936000335693359375000+67847907642187500*theta[2]... 

> nops(B); 1
 

2 

> factor(B[1]); 1
 

(8*theta[2]-25)*(7*theta[2]-500)*(theta[2]+15)*(2*theta[2]+525)*(theta[2]+75)*(theta[2]-425)*(134456*theta[2]^3-10852275*theta[2]^2-4304728125*theta[2]+935718750)
(8*theta[2]-25)*(7*theta[2]-500)*(theta[2]+15)*(2*theta[2]+525)*(theta[2]+75)*(theta[2]-425)*(134456*theta[2]^3-10852275*theta[2]^2-4304728125*theta[2]+935718750) 

> rts := [fsolve(B[1], theta[2], complex)]; 1
 

[-262.5000000, -143.2006009, -75.00000000, -15.00000000, .2172513326, 3.125000000, 71.42857143, 223.6958131, 425.]
[-262.5000000, -143.2006009, -75.00000000, -15.00000000, .2172513326, 3.125000000, 71.42857143, 223.6958131, 425.] 

Of these, note that only one is in the range  0 <= r and r <= 1. 

> fsolve(subs(theta[2] = rts[5], B[2]), theta[1]); 1
 

.5191263946 

We have exactly one root in the "probability simplex" 

 

(theta[1], theta[2]) = (.5191263946, .2172513326) 

 

Since theta[1]is relatively large here, the data would indicate 

that this region is a CG-rich area.   

>