The revisions let you track differences between multiple versions of a post.

This howto provides step-by-step instructions on writing your own objective function for OpenMx. You will first want to prototype your objective function using either an MxAlgebraObjective or an MxRObjective. An algebra objective consists of an arbitrary number of matrix algebra statements that ultimately evaluate to a 1 x 1 matrix. An R objective consists of an arbitrary R function that evaluates a model to return a numeric value. Algebra objectives are computationally faster than R objectives, but algebra objectives are more limited in their expressive power.

When you are satisfied prototyping your new objective function, we will guide you through the process of writing the front-end for this objective function in R, and writing the back-end for this objective function in C.

Can you express your objective function using the matrix operators and functions provided by OpenMx? If so, then create an MxAlgebraObjective. The following example is from the SimpleAlgebra demo.

A <- mxMatrix("Full", values = c(1, 2), nrow=1, ncol=2, name="A") B <- mxMatrix("Full", values = c(3, 4), nrow = 2, ncol = 1, free = c(TRUE, FALSE), name="B") algebra <- mxAlgebra(A %*% B - 3, "algebra") algSquared <- mxAlgebra(algebra %*% algebra, "algebraSquared") model <- mxModel("themodel") model <- mxModel(model, A, B) model <- mxModel(model, algebra, algSquared) model <- mxModel(model, mxAlgebraObjective("algebraSquared")) # Test 1: Algebra is a multiply and subtract modelOut <- mxRun(model)

Otherwise, express your objective function in the R front-end. Write a function in R that accepts two arguments. The first argument is the mxModel that should be evaluated, and the second argument is some persistent state information that can be stored between one iteration of optimization to the next iteration. It is valid for the function to simply ignore the second argument. The function must return either a single numeric value, or a list of exactly two elements. If the function returns a list, the first argument must be a single numeric value and the second element will be the new persistent state information to be passed into this function at the next iteration. The single numeric value will be used by the optimizer to perform optimization. For more help on the slots in the MxModel class, see the help provided with the OpenMx package (type '?MxModel' in R).

The following example is from the SimpleRObjective demo.

A <- mxMatrix(nrow = 2, ncol = 2, values = c(1:4), free = TRUE, name = 'A') squared <- function(x) { x ^ 2 } objFunction <- function(model, state) { values <- mxEval(A, model) return(squared(values[1,1] - 4) + squared(values[1,2] - 3) + squared(values[2,1] - 2) + squared(values[2,2] - 1)) } objective <- mxRObjective(objFunction) model <- mxModel('model', A, objective) modelOut <- mxRun(model)

In order to write your own objective function, you will need to create your own objective function class in R. For teaching purposes, we are going to use the MxAlgebraObjective class as a guide. First, your class needs to be a sub-class of the MxBaseObjective class. The MxBaseObjective class contains three slots: 'name' which is a character string, 'data' which is either a character string or a number, and 'result' which is a 1x1 matrix that stores the result of optimization. Your class should contain any additional slots that you need to pass along to the optimizer. In the MxAlgebraObjective, we have an additional slot that will inform the optimizer as to which MxAlgebra to use for optimization.

setClass(Class = "MxAlgebraObjective", representation = representation( algebra = "MxCharOrNumber"), contains = "MxBaseObjective")

Next, you will need to overwrite three generic methods for your new objective function class. These three methods are initialize, omxObjFunConvert, and omxObjFunNamespace.

The initialize method is a constructor method. It will never be directly called by the user, but it will be critical for making new instances of your class. When writing your initialize method, it is required that the default value for the data slot is as.integer(NA) and the default value for the name slot is 'objective'. In fact, the data slot will be modified only be the omxObjFunConvert method, and the name slot will only be modified by the omxObjFunNamespace method.

setMethod("initialize", "MxAlgebraObjective", function(.Object, algebra, data = as.integer(NA), name = 'objective') { .Object@name <- name .Object@algebra <- algebra .Object@data <- data return(.Object) } )

The omxObjFunConvert method is invoked to transform an objective function from a form that interacts with the R front-end, to a form that is suitable for processing by the C back-end. Any character strings in the objective function are converted into matrix/algebra pointers for the back-end to process. You should read the documentation for the function omxLocateIndex to familiarise yourself with matrix/algebra pointers. Two arguments are passed into this function, a flattened version of the model, and the original model itself. When at all possible, it is recommended that you interact with the flattened version of the model. Below we have included the omxObjFunConvert method for both the algebra objective and the FIML objective. The FIML objective has been included so that you may see how the data slot is processed.

NOTE: .Object@algebra is a string on input, and is converted into a matrix/algebra pointer.

setMethod("omxObjFunConvert", signature("MxAlgebraObjective"), function(.Object, flatModel, model) { name <- .Object@name algebra <- .Object@algebra algebraIndex <- omxLocateIndex(flatModel, algebra, name) if (is.na(algebraIndex)) { stop(paste("Could not find a matrix/algebra with name", algebra, "in the model.")) } .Object@algebra <- algebraIndex return(.Object) })

NOTE: .Object@data, .Object@means, and .Object@covariance are all strings at the beginning of this function, and are converted into data pointers or matrix/algebra pointers.

setMethod("omxObjFunConvert", signature("MxFIMLObjective"), function(.Object, flatModel, model) { name <- .Object@name if(is.na(.Object@data)) { msg <- paste("The MxFIMLObjective", omxQuotes(name), "does not have a dataset associated with it in model", omxQuotes(flatModel@name)) stop(msg, call.=FALSE) } if(flatModel@datasets[[.Object@data]]@type != 'raw') { msg <- paste("The dataset associated with MxFIMLObjective", omxQuotes(name), "in model", omxQuotes(flatModel@name), "is not raw data.") stop(msg, call.=FALSE) } if(!is.na(.Object@means)) { .Object@means <- omxLocateIndex(flatModel, .Object@means, name) } .Object@covariance <- omxLocateIndex(flatModel, .Object@covariance, name) .Object@data <- omxLocateIndex(flatModel, .Object@data, name) .Object@definitionVars <- generateDefinitionList(flatModel) return(.Object) })

The omxObjFunNamespace method is used to transform an objective function from its local model namespace, into the flattened model namespace. Basically you must use the function omxConvertIdentifier on any slot that is storing name information. Remember to always invoke omxConvertIdentifier on the 'name' slot and the 'data' slot. Shown below is the implementation for the FIML objective function.

setMethod("omxObjFunNamespace", signature("MxFIMLObjective"), function(.Object, modelname, namespace) { .Object@name <- omxIdentifier(modelname, .Object@name) .Object@covariance <- omxConvertIdentifier(.Object@covariance, modelname, namespace) .Object@means <- omxConvertIdentifier(.Object@means, modelname, namespace) .Object@data <- omxConvertIdentifier(.Object@data, modelname, namespace) return(.Object) })

Finally, you should write your own method that creates instances of your new objective function class. This interface will invoke the 'initialize' constructor that you wrote earlier, but should also include some basic sanity-checking before invoking the constructor. Here is the algebra objective function interface. Since this method does not have access to any parts of the model, you will not be able to perform full error-checking at this stage. Full error-checking should be performed in the omxObjFunConvert method.

mxAlgebraObjective <- function(algebra) { if (missing(algebra) || typeof(algebra) != "character") { stop("Algebra argument is not a string (the name of the algebra)") } return(new("MxAlgebraObjective", algebra)) }

TODO: write this section

- Login or register to post comments
- Printer-friendly version
- Login or register to post comments
- Send to friend