Good point. Now this returns 0.04062184. Hmmm..... On Mon, Oct 2, 2017 at 6:30 PM, Hollie Johnson (PGR) < h.a.john...@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-boun...@r-project.org> on behalf of Eric Berger < > ericjber...@gmail.com> > *Sent:* 02 October 2017 16:16:13 > *To:* Bert Gunter > *Cc:* r-help@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.4...@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.john...@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@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@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@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@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.
