# Issues with 'Miwa' algorithm in mvtnorm package

6 messages
Open this post in threaded view
|

## Issues with 'Miwa' algorithm in mvtnorm package

 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]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: Issues with 'Miwa' algorithm in mvtnorm package

 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) < [hidden email]> 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]] > > ______________________________________________ > [hidden email] 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]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|

## Re: 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 <[hidden email]> 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) < > [hidden email]> 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]] > > > > ______________________________________________ > > [hidden email] 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]] > > ______________________________________________ > [hidden email] 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]] ______________________________________________ [hidden email] mailing list -- To UNSUBSCRIBE and more, see https://stat.ethz.ch/mailman/listinfo/r-helpPLEASE do read the posting guide http://www.R-project.org/posting-guide.htmland provide commented, minimal, self-contained, reproducible code.
Open this post in threaded view
|