Wed, 08/24/2011 - 12:56 — neale

### Jump to:

Project: | OpenMx |

Component: | Code |

Category: | bug report |

Priority: | normal |

Assigned: | tbrick |

Status: | active |

Description

Per Hao's post http://openmx.psyc.virginia.edu/thread/1061 there is definitely odd behavior with omxMnor() with singular and non-positive definite covariance matrices. I seem to remember that Genz's integration routines were extended to handle positive semi-definite matrices (one or more eigenvalues of zero) by reducing the dimensionality of the problem. However, the answer is clearly wrong in this case:

> omxMnor(array(1,dim=c(2,2)),cbind(0,0),cbind(-Inf,-Inf),cbind(0,0))

[,1]

[1,] 0.375

(should be .5)

The behavior for negative definite matrices is less clear - we should probably throw an error, or at least flag to the optimizer that we are in infeasible region. At present it gives a value which is not correct:

> omxMnor(array(c(1,2,2,1),dim=c(2,2)),cbind(0,0),cbind(-Inf,-Inf),cbind(0,0))

[,1]

[1,] 0.4262082

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

## Comments

## #1

Swapping out SADMVN for MVNDST as the integration routine fixes the singular problem, but only pretends that the negative definite issue is also singular. No error is returned by MVNDST for this problem; we could catch it on the way in (if any correlation is >1 or <-1 then a problem).

> #allOrdinalResults <- polychoricMatrix(ordinalData,useDeviations=T)

> #summary(allOrdinalResults$estimatedModel)

> omxMnor(array(1,dim=c(2,2)),cbind(0,0),cbind(-Inf,-Inf),cbind(0,0))

[,1]

[1,] 0.5

>

>

> S<-diag(1,2);

> S[1,2]<-S[2,1]<-2

> omxMnor(S,cbind(0,0),cbind(-Inf,-Inf),cbind(0,0))

[,1]

[1,] 0.5

>

>

## #2

Unfortunately, MVNDST is MUCH slower (for the same required precision) and not really practical in an optimization context. Therefore I think we should trap non-pd before calling the integration routine.