cmlmtProfile#

Purpose#

Computes profile likelihood traces and profile t traces for models estimated using maximum likelihood.

Format#

out = cmlmtProfile(&modelProc, par[, ..., c1])#
Parameters:
  • &modelProc (function pointer) – Pointer to a procedure that computes the function to be minimized.

  • par (struct) – Instance of a PV structure, constructed using the “pack” functions.

  • ... (various) – Optional input arguments. They can be any set of structures, matrices, arrays, strings required to compute the function. Can include GAUSS data types or a DS structure for dataset manipulation. Specific usage depends on the requirements of the modelProc.

  • c1 (struct) –

    Optional input. Instance of a cmlmtControl structure containing the following members:

    c1.A

    MxK matrix, linear equality constraint coefficients: c1.A * p = c1.B where p is a vector of the parameters.

    c1.B

    Mx1 vector, linear equality constraint constants: c1.A * p = c1.B where p is a vector of the parameters.

    c1.C

    MxK matrix, linear inequality constraint coefficients: c1.C * p >= c1.D where p is a vector of the parameters.

    c1.D

    Mx1 vector, linear inequality constraint constants: c1.C * p >= c1.D where p is a vector of the parameters.

    c1.eqProc

    Scalar, pointer to a procedure that computes the nonlinear equality constraints. It has two input arguments: an instance of a PV parameter structure, and an instance of a DS data structure; and one output argument, a vector of computed equality constraints. Default = {.}, i.e., no equality procedure.

    c1.IneqProc

    Scalar, pointer to a procedure that computes the nonlinear inequality constraints. It has two input arguments: an instance of a PV parameter structure, and an instance of a DS data structure; and one output argument, a vector of computed inequality constraints. Default = {.}, i.e., no inequality procedure.

    c1.eqJacobian

    Scalar, pointer to a procedure that computes the Jacobian of the equality constraints. It has two input arguments: an instance of a PV parameter structure, and an instance of a DS data structure; and one output argument, a matrix of derivatives of the equality constraints with respect to the parameters. Default = {.}, i.e., no equality Jacobian procedure.

    c1.ineqJacobian

    Scalar, pointer to a procedure that computes the Jacobian of the inequality constraints. It has two input arguments: an instance of a PV parameter structure, and an instance of a DS data structure; and one output argument, a matrix of derivatives of the inequality constraints with respect to the parameters. Default = {.}, i.e., no inequality Jacobian procedure.

    c1.Bounds

    1x2 or Kx2 matrix, bounds on parameters. If 1x2, all parameters have the same bounds. Default = {-1e256, 1e256}.

    c1.bayesAlpha

    Exponent of the Dirichlet random variates used in the weighted bootstrap. Default = 1.4.

    c1.priorProc

    Pointer to a procedure for computing the prior. Assumes a uniform prior if not provided.

    c1.numSamples

    Number of re-samples in the weighted likelihood bootstrap.

    c1.BayesFname

    Filename for the simulated posterior parameters dataset. Defaults to a unique “BAYESxxxx” pattern.

    c1.maxBootTime

    Maximum time allowed for resampling.

    c1.algorithm

    Scalar, descent algorithm. 0 - Modified BFGS (Default), 1 = BFGS, 2 = DFP, 3 = Newton, 4 = BHHH.

    c1.useThreads

    Scalar, if nonzero, threading is turned on; otherwise, off. Default = off.

    c1.Switch

    4x1 or 4x2 vector, controls algorithm switching. If 4x1, the details follow specific conditions. If 4x2, CMLMT switches between the algorithms in column 1 and column 2. Default = {1 3, .0001 .0001, 10 10, .0001 .0001}.

    c1.lineSearch

    Scalar, sets the line search method. 0 = augmented Lagrangian penalty method (requires constraints), 1 = STEPBT (quadratic and cubic curve fit) (default), 2 = Brent’s method, 3 = BHHHStep, 4 = half, 5 = Strong Wolfe’s condition.

    c1.trustRadius

    Scalar, radius of the trust region. If missing, trust region not applied. Sets a maximum amount of the direction at each iteration. Default = .001.

    c1.penalty

    Scalar, augmentation constant for augmented Lagrangian penalty line search method.

    c1.active

    Kx1 vector, set the K-th element to zero to fix it to start value. Use the GAUSS function pvGetIndex() to determine where parameters in the PV structure are in the vector of parameters. Default = all parameters are active.

    c1.numObs

    Scalar, number of observations, required if the log-likelihood

    c1.maxIters

    Scalar, maximum number of iterations. Default = 10000.

    c1.dirTol

    Scalar, convergence tolerance. Iterations cease when all elements of the direction vector are less than this value.

    c1.weights

    Vector, weights for objective function returning a vector. Default = 1.

    c1.CovParType

    Scalar, if 2, QML covariance matrix; else if 0, no covariance matrix is computed; else ML covariance matrix is computed. Default = 1.

    c1.alpha

    Scalar, probability level for statistical tests. Default = .05.

    c1.FeasibleTest

    Scalar, if nonzero, parameters are tested for feasibility before computing function in line search. If function is defined outside inequality boundaries, then this test can be turned off. Default = 1.

    c1.MaxTries

    Scalar, maximum number of attempts in random search. Default = 100.

    c1.randRadius

    Scalar, if zero, no random search is attempted. If nonzero, it is the radius of the random search. Default = .001.

    c1.gradMethod

    Scalar, method for computing numerical gradient. 0 = central difference, 1 = forward difference (default), 2 = backward difference, 3 = complex derivatives.

    c1.hessMethod

    Scalar, method for computing numerical Hessian. 0 = central difference, 1 = forward difference (default), 2 = backward difference.

    c1.gradStep

    Scalar or Kx1, increment size for computing numerical gradient. If scalar, stepsize will be value times parameter estimates for the numerical gradient. If Kx1, the step size for the gradient will be the elements of the vector.

    c1.hessStep

    Scalar or Kx1, increment size for computing numerical Hessian. If scalar, stepsize will be value times parameter estimates for the numerical Hessian. If Kx1, the step size for the Hessian will be the elements of the vector.

    c1.gradCheck

    Scalar, if nonzero and if analytical gradients and/or Hessian have been provided, numerical gradients and/or Hessian are computed and compared against the analytical versions.

    c1.state

    Scalar, seed for random number generator.

    c1.title

    String, title of run.

    c1.PrintIters

    Scalar, if nonzero, prints iteration information. Default = 0.

    c1.disableKey

    Scalar, if nonzero, keyboard input disabled.

Returns:

out (struct) –

An instance of a cmlmtResults structure. Contains the results of the optimization problem, including parameter estimates, function evaluations, and various statistical measures.

out.par

Instance of a PV structure containing the parameter estimates will be placed in the member matrix out.par.

out.fct

Scalar, function evaluated at parameters in par.

out.returnDescription

String, description of return values.

out.covPar

KxK matrix, covariance matrix of parameters.

out.covParDescription

String, description of covPar.

out.numObs

Scalar, number of observations.

out.hessian

KxK matrix, Hessian evaluated at parameters in par.

out.xproduct

KxK matrix, cross-product of NxK matrix of first derivatives evaluated at parameters in par. Not available if log-likelihood function returns a scalar.

out.waldLimits

Kx2 matrix, Wald confidence limits.

out.inverseWaldLimits

Kx2 matrix, confidence limits by inversion of Wald statistics. Available only if cmlmtInverseWaldLimits() has been called.

out.profileLimits

Kx2 matrix, profile likelihood confidence limits, i.e., by inversion of likelihood ratio statistics. Only available if cmlmtProfileLimits() has been called.

out.bayesLimits

Kx2 matrix, weighted likelihood Bayesian confidence limits. Only available if cmlmtBayes() has been called.

out.bootLimits

Kx2 Matrix, bootstrap confidence limits. Available only if cmlmtBoot() has been called.

out.gradient

Kx1 vector, gradient evaluated at the parameters in par.

out.numIterations

Scalar, number of iterations.

out.elapsedTime

Scalar, elapsed time of iterations.

out.alpha

Scalar, probability level of confidence limits. Default = .05.

out.title

String, title of run.

out.lagr

An instance of a cmlmtLagrange structure containing the Lagrangeans for the constraints. For an instance named lagr, the members are:

  • out.lagr.lineq - Mx1 vector, Lagrangeans of linear equality constraints,

  • out.lagr.nlineq - Nx1 vector, Lagrangeans of nonlinear equality constraints

  • out.lagr.linineq - Px1 vector, Lagrangeans of linear inequality constraints

  • out.lagr.nlinineq - Qx1 vector, Lagrangeans of nonlinear inequality constraints

  • out.lagr.bounds - Kx2 matrix, Lagrangeans of bounds

  • out.Lagr.EqCov - (M+N)x(M+N) matrix, covariance matrix of equality constraints

  • out.Lagr.IneqCov - (P+Q)x(P+Q) matrix, covariance matrix of inequality constraints

Whenever a constraint is active, its associated Lagrangean will be nonzero. For any constraint that is inactive throughout the iterations as well as at convergence, the corresponding Lagrangean matrix will be set to a scalar missing value.

out.retcode

Return code:

  • 0 - Normal convergence

  • 1 - Forced exit

  • 2 - Maximum number of iterations exceeded

  • 3 - Function calculation failed

  • 4 - Gradient calculation failed

  • 5 - Hessian calculation failed

  • 6 - Line search failed

  • 7 - Functional evaluation failed

  • 8 - Error with initial gradient

  • 9 - Error with constraints

  • 10 - Second update failed

  • 11 - Maximum time exceeded

  • 12 - Error with weights

  • 13 - Quadratic program failed

  • 14 - Equality constraint Jacobian failed

  • 15 - Inequality constraint Jacobian failed

  • 16 - Function evaluated as complex

  • 20 - Hessian failed to invert

  • 34 - Data set could not be opened

Example#

Biochemical Oxygen Demand (BOD) Analysis#

This example demonstrates the application of cmlmtProfile() to model the Biochemical Oxygen Demand (BOD) data using a log-likelihood function.

library cmlmt;

// Define the log-likelihood function
proc lnlk(struct PV p, struct DS d, ind);
    local dev, s2, m, r, b0, b;

    // Declare 'mm' to be a modelResults
    // struct local to this procedure, 'lnlk'
    struct modelResults mm;

    // Unpack parameters
    b0 = pvUnpack(p, 1);
    b = pvUnpack(p, 2);

    // Calculate model predictions
    r = exp(-b * d[2].dataMatrix);
    m = 1 - r;

    // Calculate deviations
    dev = d[1].dataMatrix - b0 * m;
    s2 = dev'dev/rows(dev);

    // Calculate log-likelihood
    if ind[1];
        mm.function = lnpdfmvn(dev, s2);
    endif;

    // Calculate gradient, if requested
    if ind[2];
        mm.gradient = (dev / s2) .* (m ~ b0 * d[2].dataMatrix .* r);
    endif;

    retp(mm);
endp;

// Data setup
struct DS d0;
d0 = reshape(dsCreate, 2, 1);
d0[1].dataMatrix = {8.3, 10.3, 19.0, 16.0, 15.6, 19.8};
d0[2].dataMatrix = {1, 2, 3, 4, 5, 7};

// Parameter setup
struct PV p0;
p0 = pvPacki(pvCreate, 19.143, "b0", 1);
p0 = pvPacki(p0, .5311, "b", 2);

// Control structure setup
struct cmlmtControl c0;
c0 = cmlmtControlCreate;
c0.Bounds = {10 35, 0 2};  // Set parameter bounds

// Perform the profile likelihood analysis
struct cmlmtResults out;
out = cmlmtProfile(&lnlk, p0, d0, c0);

Remarks#

  • cmlmtProfile() is utilized to explore the parameter space of maximum likelihood estimates more thoroughly, offering insights into the confidence intervals and sensitivity of the estimates.

  • This function is especially useful in complex models where the standard error may not provide a complete picture of parameter uncertainty.

  • The control structure allows extensive customization of the profiling process, making it adaptable to a wide range of models and research questions.