Suppose I am running a Monte Carlo, and in each simulation, CMLMT is called. Most of the time, CMLML works find, but there are occasions when CML may fail, for example, due to "Not implemented for complex matrices." If this happens, I would like to discard this simulation and start a new simulation. Typically, I do not know what type of error CMLMT may create, as the simulated data may be ill-behaved in many different ways. Is there a decent way of "trapping" such that whenever CMLMT generates errors (and the simulation will stop if without "trapping"), the simulation loop will skip this round of simulation and continue to the next simulation?

## 3 Answers

0

If these errors are happening inside of your likelihood function, you could check these inputs to see if they are complex or Nans. For example, if we change the likelihood procedure in `cmlmt1.e` to artificially create a complex number before it is handed off to `cdfn` which does not accept complex numbers like this:

proc lpr(struct PV p, struct DS d, ind); local s2,b0,b,y,x,yh,u,res,g1,g2; struct modelResults mm; b0 = pvUnpack(p,"b0"); b = pvUnpack(p,"b"); s2 = pvUnpack(p,"variance"); y = d[1].dataMatrix; x = d[2].dataMatrix; yh = b0 + x * b; res = y - yh; u = y[.,1] ./= 0; if ind[1]; mm.function = u.*lnpdfmvn(res,s2) + (1-u).*(ln(cdfnc(yh/sqrt(s2)))); endif; if ind[2]; //make negative so 'yh' will be complex s2 = -s2; yh = yh/sqrt(s2); //stop program if 'yh' is complex if iscplx(yh); stop; endif; g1 = ((res~x.*res)/s2)~((res.*res/s2)-1)/(2*s2); g2 = ( -( ones(rows(x),1)~x )/sqrt(s2) )~(yh/(2*s2)); g2 = (pdfn(yh)./cdfnc(yh)).*g2; mm.gradient = u.*g1 + (1-u).*g2; endif; retp(mm); endp;

Now with this new likelihood procedure, the example file `cmlmt1.e` will catch the problem before it causes an error and the program will `stop`. Once we have done this, there is still one problem. You don't want the entire program to stop.

You can accomplish this by running this file with the `run` command using the `-r` flag. For example, if you made a main file which we will call `cmlmtdriver.gss` with these contents:

for i(1, 3, 1); print "starting iteration " i; run -r cmlmt1.e; endfor; print "main file still running";

If you have any questions about how to adapt these techniques to your simulation program, please ask.

`Other useful functions`

Besides just checking whether or not a matrix is complex with `iscplx`, you can check whether or not a matrix has a complex portion that is larger than a certain tolerance with the `hasimag` function. This will allow you to retain only the real portion of the matrix if the complex portion is 0 to machine precision.

The following GAUSS functions may also be helpful in this context: `ismiss`, `isinfnanmiss`, `scalmiss`.

There is also the `trap` function which can be used with `chol`, `invpd`, or `solpd` to trap in the case that the input matrix is not positive definite. `trap` can also be used with the slash operator and `inv` to trap cases in which the input matrix is singular, or with the eigenvalue functions to trap the case in which the eigenvalues cannot be calculated.

0

Thanks!

My code indicates that "cholup" and "choldn" sometimes give "error G0043 : Not implemented for complex matrices." The code from CMLMT.src has the following:

` oldt = trapchk(1);`

trap 1,1;

...

code involving cholup/choldn

...

trap oldt,1;

It appears that error trapping does not here for "cholup" and "choldn."

0

You are correct, `cholup` and `choldn` do not support the `trap` command.

One thing to know is that if a complex matrix is sent to `cholup` or `choldn`, they do perform a check on the complex portion of the input to see if it is significant like the `hasimag` function.

It looks like all the calls to `cholup` and `choldn` in CMLMT take the hessian as an input. I think the first step I would take is to add an `iscplx` check to any outputs from your likelihood function and add a `stop` if the output was complex.

If the programs take a while to run, I might also add an `if iscplx()` check to some of the other variables being passed into `choldn` or `cholup`, particularly `g0` and `v2`. Or if you just want to be done with it, you could take a line like this:

h1 = cholup(h1,g1/v3-g0/v3+ak*dx/v3);

and change it to something like this:

local cholup_input; cholup_input = g1/v3-g0/v3+ak*dx/v3; if iscplx(h1) or iscplx(cholup_input); errorlogat "exiting CMLMT due to unexpected complex matrix"; stop; endif; h1 = cholup(h1,cholup_input);

## Your Answer

## 3 Answers

If these errors are happening inside of your likelihood function, you could check these inputs to see if they are complex or Nans. For example, if we change the likelihood procedure in `cmlmt1.e` to artificially create a complex number before it is handed off to `cdfn` which does not accept complex numbers like this:

proc lpr(struct PV p, struct DS d, ind); local s2,b0,b,y,x,yh,u,res,g1,g2; struct modelResults mm; b0 = pvUnpack(p,"b0"); b = pvUnpack(p,"b"); s2 = pvUnpack(p,"variance"); y = d[1].dataMatrix; x = d[2].dataMatrix; yh = b0 + x * b; res = y - yh; u = y[.,1] ./= 0; if ind[1]; mm.function = u.*lnpdfmvn(res,s2) + (1-u).*(ln(cdfnc(yh/sqrt(s2)))); endif; if ind[2]; //make negative so 'yh' will be complex s2 = -s2; yh = yh/sqrt(s2); //stop program if 'yh' is complex if iscplx(yh); stop; endif; g1 = ((res~x.*res)/s2)~((res.*res/s2)-1)/(2*s2); g2 = ( -( ones(rows(x),1)~x )/sqrt(s2) )~(yh/(2*s2)); g2 = (pdfn(yh)./cdfnc(yh)).*g2; mm.gradient = u.*g1 + (1-u).*g2; endif; retp(mm); endp;

Now with this new likelihood procedure, the example file `cmlmt1.e` will catch the problem before it causes an error and the program will `stop`. Once we have done this, there is still one problem. You don't want the entire program to stop.

You can accomplish this by running this file with the `run` command using the `-r` flag. For example, if you made a main file which we will call `cmlmtdriver.gss` with these contents:

for i(1, 3, 1); print "starting iteration " i; run -r cmlmt1.e; endfor; print "main file still running";

If you have any questions about how to adapt these techniques to your simulation program, please ask.

`Other useful functions`

Besides just checking whether or not a matrix is complex with `iscplx`, you can check whether or not a matrix has a complex portion that is larger than a certain tolerance with the `hasimag` function. This will allow you to retain only the real portion of the matrix if the complex portion is 0 to machine precision.

The following GAUSS functions may also be helpful in this context: `ismiss`, `isinfnanmiss`, `scalmiss`.

There is also the `trap` function which can be used with `chol`, `invpd`, or `solpd` to trap in the case that the input matrix is not positive definite. `trap` can also be used with the slash operator and `inv` to trap cases in which the input matrix is singular, or with the eigenvalue functions to trap the case in which the eigenvalues cannot be calculated.

Thanks!

My code indicates that "cholup" and "choldn" sometimes give "error G0043 : Not implemented for complex matrices." The code from CMLMT.src has the following:

` oldt = trapchk(1);`

trap 1,1;

...

code involving cholup/choldn

...

trap oldt,1;

It appears that error trapping does not here for "cholup" and "choldn."

You are correct, `cholup` and `choldn` do not support the `trap` command.

One thing to know is that if a complex matrix is sent to `cholup` or `choldn`, they do perform a check on the complex portion of the input to see if it is significant like the `hasimag` function.

It looks like all the calls to `cholup` and `choldn` in CMLMT take the hessian as an input. I think the first step I would take is to add an `iscplx` check to any outputs from your likelihood function and add a `stop` if the output was complex.

If the programs take a while to run, I might also add an `if iscplx()` check to some of the other variables being passed into `choldn` or `cholup`, particularly `g0` and `v2`. Or if you just want to be done with it, you could take a line like this:

h1 = cholup(h1,g1/v3-g0/v3+ak*dx/v3);

and change it to something like this:

local cholup_input; cholup_input = g1/v3-g0/v3+ak*dx/v3; if iscplx(h1) or iscplx(cholup_input); errorlogat "exiting CMLMT due to unexpected complex matrix"; stop; endif; h1 = cholup(h1,cholup_input);