[R] Issues with 'Miwa' algorithm in mvtnorm package

[Bert Gunter]

Rather specialized.

As this appears to be primarily a statistical, not an R programming
question, you may do better posting on a statistical site like
http://stats.stackexchange.com/ if you don't get a satisfactory reply here
. Alternativey, if you think this is a package bug, perhaps contact the
package maintainer directly, as (s)he may not monitor this list.

-- Bert



Bert Gunter

"The trouble with having an open mind is that people keep coming along and
sticking things into it."
-- Opus (aka Berkeley Breathed in his "Bloom County" comic strip )

On Mon, Oct 2, 2017 at 5:16 AM, Hollie Johnson (PGR) <
h.a.johnson at newcastle.ac.uk> wrote:

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

	[[alternative HTML version deleted]]