Marginal Likelihood Integrals for Mixtures of Independence Models

Marginal Likelihood Integrals
for Mixtures of Independence Models

Shaowei Lin, Bernd Sturmfels, and Zhiqiang Xu



View this page in Romanian courtesy of azoft
View this page in Polish courtesy of Valeria Aleksandrova
View this page in Danish courtesy of ScienceLakes

Research Paper

Our paper "Marginal Likelihood Integrals for Mixtures of Independence Models" can be found here.

Maple Library

The Maple library can be downloaded from here. To use it, put the file in your work directory and type the command
> read "integrals.lib";
when you start your maple session.

Toric Matrix

Suppose our independence model is defined by parameters s1, ..., sk and t1, ..., tk. Let s and t be lists storing the parameters. The function
ToricMatrix(s,t)
will produce an array storing the reduced matrix A, i.e. repeated columns are removed. Note that the columns of A are sorted lexicographically. For instance, the matrix A in Example 2.1 is given by
> A:=ToricMatrix([1,2],[1,1]);

[1    1    1    0    0    0]
[                          ]
[0    0    0    1    1    1]
[                          ]
[2    1    0    2    1    0]
[                          ]
[0    1    2    0    1    2]
To get the non-reduced matrix A, use the option "reduced=false". The following matrix is from Example 2.5.
> A:=ToricMatrix([4],[1],reduced=false);

[4    3    3    3    3    2    2    2    2    2    2    1    1    1    1    0]
[                                                                            ]
[0    1    1    1    1    2    2    2    2    2    2    3    3    3    3    4]

Marginal Likelihood for Mixture Model

Let U be a list storing the data vector. Make sure that the i-th entry of U corresponds with the i-th column of the matrix A produced by the function ToricMatrix. Then, the marginal likelihood for the data vector U with respect to the mixture model is given by
ML(s,t,U)
where the output depends on whether U is a reduced or non-reduced data vector. To understand the difference between the two, we study the coin toss example in Example 2.5. Suppose we have the data vector
(U0000,U0001,U0010,...,U1111) = (51,4,4,5,5,12,12,12,12,12,13,6,6,6,7,75).
Note that the reduced form of this data vector is
(U0,U1,U2,U3,U4) = (51,18,73,25,75).
Then, the marginal likelihood of the non-reduced data vector (U0000,U0001,...,U1111) is
ML([4],[1],[51,4,4,5,5,12,12,12,12,12,13,6,6,6,7,75])
which is approximately 0.2857408184x10-30. It is given by the formula
                   242!
---------------------------------------- I
51!4!4!5!5!12!12!12!12!12!13!6!6!6!7!75!
where I is the integral (21) in the paper. Meanwhile, the marginal likelihood of the reduced data vector (U0,U1,U2,U3,U4) is the sum of all marginal likelihoods of data vectors (U0000,U0001,...,U1111) which reduce to the given reduced data vector. We compute it with
ML([4],[1],[51,18,73,25,75])
which is approximately 0.7788716339x10-22. It is given by the formula
       242!
--------------- 418 673 425 I
51!18!73!25!75!
Note that the difference between the two marginal likelihoods is only due to the normalizing constants.

Dirichlet Priors

To compute marginal likelihoods with Dirichlet priors α, β and γ, we use the function
MLDir(s,t,U,a,b,g)
where a, b and g are lists representing α, β and &gamma respectively. Note that a, b and g are of length 2, d and d respectively where d is the number of rows in the toric matrix A, and that their entries must be positive.

Marginal Likelihood for Independence Model

The marginal likelihood for the independence model can be computed with the option "mixed=false".
ML(s,t,U,mixed=false)
MLDir(s,t,U,[],b,[],mixed=false)
Note that for the function MLDir, one needs to put dummy lists such as "[]" in the place of the the parameters a and g, as shown above.

Bayes Factor

To compute the Bayes factor mentioned in Corollary 5.4 of the paper, simply evaluate
ML(s,t,U)/ML(s,t,U,mixed=false)
MLDir(s,t,U,a,b,g)/MLDir(s,t,U,[],b,[],mixed=false)

Other Options

The functions ML and MLDir have the following additional options:
  1. symbolic = true / false (default: false)
    Tells the function whether to use the symbolic expansion algorithm described in Section 4.2 of the paper. This algorithm may be faster when the expanded integrand has relatively few terms.
  2. float = true / false (default: false)
    Tells the function whether to compute the integral with floating point arithmetic or with exact rational number arithmetic. Use the variable Digits to control the precision of the calculations. For instance, the command Digits := 50; sets the precision to 50 decimal digits. Note that when computing marginal likelihoods with non-integer Dirichlet prior parameters, the float option is automatically set to true.
  3. quiet = true / false (default: true)
    Tells the function to withhold or print out the progress of the computation respectively. Useful when the computation takes several hours or days. Use the command interface(quiet = true) to turn off the system's memory and time print-out.
  4. jump = positive integer (default: 1000)
    Used in conjunction with the option quiet = false. This variable tells the function how frequent the progress print-outs should be made. It is the number of terms of the integral formula which are summed up before a progress print-out is made.

Convenience Functions

The library has some useful convenience functions.
  1. ToricNumRows(s,t), ToricNumCols(s,t)
    Returns d and n, the number of rows and cols respectively for the reduced matrix of the independence model defined by parameters s and t. Use the option "reduced=false" to compute the same for the non-reduced matrix.
  2. ReducedVector(s,t,U)
    Returns the reduced data vector for an unreduced data vector U of the independence model with parameters s and t.
  3. NumTerms(A,U)
    Returns the number of terms in the expansion of the integrand defined by A and U.
  4. NumTermsBound(A,U)
    Returns a list [l,u] where l and u are the lower and upper bounds given in Theorem 3.5 for the number of terms in the expansion of the integrand defined by A and U.
  5. ExpandHalfPoly(A,U)
    Because of Proposition 3.2(1), only half the terms of every polynomial are stored. This function expands the integrand defined by A and U, and returns the result in an array as a "half-polynomial". Each row of the array corresponds to a term in the expansion, and the columns store the exponents and coefficient of each term. In particular, the coefficients are stored in the last column.
  6. GetPoly(f,[var list])
    This function takes an array f storing a half-polynomial as described previously, and returns the standard full polynomial that f represents in the variables [var list]. For example,
    > A:= ToricMatrix([2],[1]):
    > f:= ExpandHalfPoly(A,[1,1,1]):
    > GetPoly(f,[x,y]);

              2          3      3  3      3    2  2    2 
    g := 1 + x  + x y + x  y + x  y  + x y  + x  y  + y   

    > factor(g);

    2        2            
    (y  + 1) (x  + 1) (x y + 1)
  7. GetHalfPoly(f,[var list],A,U)
    This function takes an full polynomial f coming from the expansion of the integrand defined by the matrix A and data U, and returns an array storing the half polynomial. The variables in [var list] must be single-character alphabets a, b, c, d, ..., w, x, y, z.
  8. NormConst(s,t,U)
    Calculates the normalizing constant for the marginal likelihood of the data U in the mixture model with parameters s and t. In particular, it checks if U is a reduced data vector and returns the corresponding normalizing constant.
  9. ArrayNumRows(A), ArrayNumCols(A)
    Takes an array A and returns the number of rows and the number of cols it has respectively. For instance, you can use these functions to find the number of terms and variables in a half polynomial.

Examples

The following are some examples from the paper:
  1. Coin Toss (see Example 2.1)
    ML([4],[1],[51,18,73,25,75]);
  2. 100 Swiss Francs (see Section 1)
    ML([1,1],[3,3],[4,2,2,2,2,4,2,2,2,2,4,2,2,2,2,4]);
  3. Schizophrenic Patients (see Section 5)
    ML([1,1],[2,2],[43,16,3,6,11,10,9,18,16]);

Suggestions and Bugs

If you find any bugs in the code or have any comments or suggestions, please do not hesitate to contact Shaowei Lin. His email address is
shaowei [at] math [dot] berkeley [dot] edu


Last Updated: 14 Sep 2008