Tue, 12/08/2009 - 15:22

A recently posted issue requested the implementation of Mx1's `\mnor`

and `\allint`

functions, which perform integrations between thresholds on a multivariate normal distribution. This functionality is already being used in FIML Objectives if the data are ordinal (using code supplied by Alan Genz), so the back-end implementation is not difficult.

The implementation does raise some front-end questions, however. The Mx1 version of these functions use an augmented matrix structure (described below) as input. I'm interested in feedback about how it should work in OpenMx.

The Mx1 Implementation:

-------------------------------

`\mnor()`

calculates the integral between upper and lower limits on each variable (that is, the volume of a hypercube between the specified limits). For `\mnor()`

the input matrix is the vertical adhesion of:

- Covariance Matrix
- Vector of Means
- Vector of Upper Integration Limit
- Vector of Lower Integration Limit
- Vector containing a code for the type of integration required (that is, which of the limits above should be +INF or -INF)

`\mnor()`

returns a 1x1 matrix containing the integral.

`\allint()`

calculates the probabilities of all cells in a multivariate normal distribution that has been "sliced" by a varying number of hyperplanes along each dimension. For `\allint()`

, the input is:

- Covariance Matrix
- Vector of Means
- Vector Containing the Number of Thresholds/Hyperplanes
- Threshold Matrix (locations of the hyperplanes)

`\allint()`

in Mx1 returns a vector of values representing the integrals of each cell in the space, with the dimension represented by the rightmost covariance column cycling most rapidly.

Comments:

--------------

I would rather avoid the use of augmented matrices and the use of special codes if there's an easier way. I'm inclined to say that our implementation of `\mnorm()`

could simply take a matrix and vectors for the upper and lower limits, and just use the special values for infinity (or use NA) as the limits when needed.

The catch is that avoiding augmented matrices means that each argument passed in is functionally different from the others. That is, the first argument is assumed to be a covariance matrix, the second a means vector, etc. We've avoided this kind of thing in OpenMx so far, but it seems like more Rish way of doing it than augmented matrices.

`\allint`

provides a bigger problem, since there might be different numbers of thresholds for each column of the covariance matrix, so a matrix representation becomes less appropriate. I'd propose that `\allint()`

take a covariance matrix, a vector of means, and then a vector of thresholds for each column in the covariance matrix. While this makes conformance checking a little strange, it allows the user to specify any number of thresholds for a given column. What `\allint()`

outputs is also unclear, since the result should theoretically have as many dimensions as the distribution does. We'll have to flatten it, which might be ugly.

We've been keeping to the associated R functions so far--we could use the format of function `pmvnorm(lower, upper, mean, correlation/covariance)`

from the mvtnorm package to sub in for `\mnor()`

, but there's no associated function for `\allint()`

.

So, the questions:

- How would people like to specify these functions?
- Are we comfortable having functionally different arguments in OpenMx Algebras? We've avoided this so far, but it might be helpful here.
- Should we follow R-based function names for less-often-used functions, stick to Mx1 names, or branch out on our own?

Great!

I agree that augmented matrices are not very pretty. They were merely a convenient way of bundling arguments in a kludgy front end, aka Mx1.

mxIntNorm(Sigma,Mu,Lower,Upper)

works for me, with +inf and -inf as allowable arguments. Some helper function to parse an Mx1 function call would be nice, since I have scripts with many \mnorm() calls.

mxAllIntNorm(Sigma,Mu,Thresholds)

would also work, as long as the threshold matrix is padded with NA's for variables that have fewer than the maximum number of thresholds of any variable. The function could figure out how many thresholds in each dimension by finding the first NA in the corresponding column of the threshold matrix. For return, a matrix of dimension 1 x nint would be fine. For compatibility, these should be sorted so that the last variable (last column in Sigma/Mu/Threshold matrix) changes most rapidly, then the penultimate, then the penpenultimate :) etc..

One minor comment that is orthogonal to tbrick's issues: whatever functions are added to the back-end for MxAlgebra expressions, a function with the same name must be implemented in R. Otherwise conformability testing of the algebras in the front-end will throw an error.