Hollie Johnson (PGR)

2017-Oct-02 15:30 UTC

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

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