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

[Eric Berger]

Good point. Now this returns 0.04062184. Hmmm.....

On Mon, Oct 2, 2017 at 6:30 PM, Hollie Johnson (PGR) <
h.a.johnson at newcastle.ac.uk> wrote:

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