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

[Hollie Johnson (PGR)]

Hi Eric,


Thanks for having a look into this. I think you have a small typo...

B <- matrix(x, nrow=3, byrow = TRUE) should read B <- matrix(y, nrow=3, byrow = TRUE)


Regards, Hollie

________________________________
From: R-help <r-help-bounces at r-project.org> on behalf of Eric Berger <ericjberger at gmail.com>
Sent: 02 October 2017 16:16:13
To: Bert Gunter
Cc: r-help at r-project.org; Hollie Johnson (PGR)
Subject: Re: [R] Issues with 'Miwa' algorithm in mvtnorm package

Hi Hollie,
I tried to reproduce your example but could not. Here is what I did. Please
explain if you did something different.

x <- c(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)
A <- matrix(x, nrow=3, byrow = TRUE)

y <- c(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)

B <- matrix(x, nrow=3, byrow = TRUE)

pmvnorm(
  mean = rep(0, 3),
  lower = rep(-Inf, 3),
  upper = rep(0, 3),
  sigma =  A,
  algorithm = 'Miwa'
) [1]

# -10.76096       <--  this is the output from the above command

pmvnorm(
  mean = rep(0, 3),
  lower = rep(-Inf, 3),
  upper = rep(0, 3),
  sigma =  B,
  algorithm = 'Miwa'
) [1]

# -10.76096       <--  this is the output from the above command

i.e. the outputs agree in the two cases.

Regards,
Eric




On Mon, Oct 2, 2017 at 6:03 PM, Bert Gunter <bgunter.4567 at gmail.com> wrote:

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

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