You are here

omxMnor() problem

Per Hao's post https://openmx.ssri.psu.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

Reporter: 
Created: 
Wed, 08/24/2011 - 13:56
Updated: 
Thu, 12/03/2015 - 13:53

Comments

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
>
>

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.

This appears to be a duplicate issue and is generally resolved. The case of zero eigenvalues is not really correct this, but correlations greater than one are caught.