// Start by modifying the path were the library is // located in your home directory // In my case, the library is in the path LIB "~/likelihood/software/maxlike.lib"; /////////////////////////////////////////////////// // // // Parametric Maximization // // // /////////////////////////////////////////////////// // 1. // First we create the ring we will be working with ring r = 0,(x,y),dp; // Now we input the polynomial parametrization of the distribution poly f1,f2 = 7x+13y+23, -5x+7y-8; // We consider these polynomials as an ideal ideal I = f1,f2,1-f1-f2; // We define the data set as an integer vector intvec u = 11,23,31; // critical is a function that computes the critical ideal ideal J = critical(I,u); // Info returns two values (codimension, degree) Info(J); // The function MLeval has 4 parameters (the last one is optional). // The ideal I corresponding to the polynomial parametrization, // the critical ideal J, the data u, and // an optional parameter(s) corresponding to the options // for the "solve" command. In this example, the last option // stands for the precision of output in digits, which is 15. // MLeval returns a SORTED list whose entries are triplets of the form // (a,b,c) where a is a root of the critical ideal, // b is the evaluation of the distribution at the MLE a // that is b = (f1(a),f2(a),f3(a),f4(a)), // c is the value of the ML at a, i.e., c = \prod_i f_i^{u_i} // Note: MLeval changes the ambient ring to a new ring called rC MLeval(I,J,u,15); // We can check the list of all root of J with the following command setring r; Lroots(J); // We can check the list of all real roots of J setring r; Lrealroots(J); /////////////////////////////////////////////////////// /////////////////////////////////////////////////////// // We can also go directly to the MLE computation // from the ideal I and the data u, as shown in //the next example // 2. ring r = 0,(x,y),dp; poly f1,f2 = 7x+13y+23, -5x+7y-8; ideal I = f1,f2,1-f1-f2; intvec u = 11,23,31; // Since the critical ideal is 0 dimensional we can do all // the steps in the algorithm at once with the command // MLE which takes as an input the ideal corresponding to // the polynomial parametrization and the data set // (and an optional parameter(s) just as in MLeval. // The output consists of a list with two elements // the first element is equal to the output of MLeval. // The second element is the ML degree. MLE(I,u); ///////////////////////////////////////////////////////// ///////////////////////////////////////////////////////// // 3. First example in [Hosten-Khetan-Sturmfels] // The important part of this example is that // the critical ideal is not 0-dimensional // so we cannot use the command MLE directly. // First, we need to saturate the critical ideal // as explained in [Hosten-khetan-Sturmfels] ring r = 0,(p,s,t),dp; poly f1 = p*(1-s)^4 + (1-p)*(1-t)^4; poly f2 = 4*p*s*(1-s)^3 + 4*(1-p)*t*(1-t)^3; poly f3 = 6*p*s^2*(1-s)^2 + 6*(1-p)*t^2*(1-t)^2; poly f4 = 4*p*s^3*(1-s) + 4*(1-p)*t^3*(1-t); poly f5 = p*s^4 + (1-p)*t^4; ideal I = f1,f2,f3,f4,f5; intvec u = 3,5,7,11,13; ideal J = critical(I,u); J = sat(J,s-t)[1]; J = sat(J,p)[1]; J = sat(J,p-1)[1]; Info(J); list l = MLeval(I,J,u,15); l; /////////////////////////////////////////////////// // // // Constrained Maximization // // // /////////////////////////////////////////////////// // 4. First example in [Hosten-Khetan-Sturmfels] // Implicit Version // This example repeats the computation of example 3. // But this time we use the constrined optimization // approach starting from the implicit ideal instead // of the parametrization. ring r = 0,(p0,p1,p2,p3,p4),dp; matrix P[3][3] = 12*p0, 3*p1, 2*p2, 3*p1, 2*p2, 3*p3, 2*p2, 3*p3, 12*p4; ideal I = det(P); intvec u = 51,18,73,25,75; // The following function computes the Likelihood ideal // of I given the data u ideal J = Lideal(I,u); Info(J); // The previous line shows that we neet to saturate by // the singular locus since the degree of J is 13 instead // of 12 which is the ML degree as reported in [HKS] ideal Q = slocus(I); J = sat(J,Q)[1]; // Now, we need to add the condition that the sum of all // probabilities is 1 ideal m = maxideal(1); J = J,sum(m)-1; Info(J); // Finally, we use the function MLeval as explained // in the previous examples. list l = MLeval(m,J,u); l; // 5. [Example 10, HKS] ring r = 0,(p00,p01,p02,p10,p11,p12,p20,p21,p22),dp; ideal m = maxideal(1); matrix P[3][3] = m; print(P); // This time we use the command minor() instead of // det() just to show how to obtain an ideal of minors // in Singular ideal I = minor(P,3); intvec u = 16,17,7,18,3,12,1,8,16; ideal J = Lideal(I,u); Info(J); ideal Q = slocus(I); J = sat(J,Q)[1]; poly sumall = sum(m)-1; J = J+sumall; Info(J); list l = MLeval(m,J,u); l; // One can ask for all roots of the // Likelihood ideal setring r; Lroots(J); // or the real root setring r; Lrealroots(J);