Currently doing some work on local maxima on a random field and have
encountered an issue with the Miwa algorithm used with the pmvnorm function in
the mvtnorm R package.
Based on recommendations by Mi et al., we ran the mvtnorm package using the
Miwa algorithm, since we have a maximum of 4 dimensions with non-singular
matrices. However, running the estimation procedure in this way, we obtained
inconsistent results. For example, using matrices A and B, it is clear to see
that matrix B is the results of exchanging position 1 with position 3 in matrix
A.
A =
9.358*10^-3 -8.165*10^-3 -1.689*10^-8
-8.165*10^-3 9.358*10^-3 1.854*10^-8
-1.689*10^-8 1.854*10^-8 9.358*10^-3
B =
9.358*10^-3 1.854*10^-8 -1.689*10^-8
1.854*10^-8 9.358*10^-3 -8.165*10^-3
-1.689*10^-8 -8.165*10^-3 9.358*10^-3
The determinants of both matrices are small but equal, so we would expect any
issues arising from the matrix being 'close' to singular to be apparent in both
cases. The table below shows the output from pmvnorm calculated using the two
matrices A and B, and the two different algorithms, Miwa and GenzBretz using
the code below:
pmvnorm(
mean = rep(0, 3),
lower = rep(-Inf, 3),
upper = rep(0, 3),
sigma = A,
algorithm = 'Miwa'
)[1]
The results are as expected, except when using matrix A with the Miwa algorithm.
Matrix Miwa GenzBrentz
--------------------------------------
A -10.766 0.041
B 0.041 0.041
--------------------------------------
Further investigation indicates that it is the values in locations (1,3) and
(3,1) which are causing the issues; any values in the range (5*10^-7, 5*10^-9)
and (-5*10^-9, -5*10^-7) give unusual results. Can anyone help?
[[alternative HTML version deleted]]
______________________________________________
R-help@r-project.org mailing list -- To UNSUBSCRIBE and more, see
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.
