Hollie Johnson (PGR)
2017-Oct-02 12:16 UTC
[R] 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]]
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]]
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]]
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]]
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]]
Hi Hollie,
Clearly a negative number for the result of pmvnorm() makes no sense.
On the plus side, the routine seems to know it failed since the
"error"
attribute returns an NA.
You can check for this error condition programmatically, as follows:
myprob <- pmvnorm(mean=rep(0, 3),lower=rep(-Inf, 3), upper=rep(0, 3),
sigma=A, algorithm=Miwa)
if ( is.na(attr(myprob,"error"))) stop("calculation failed")
By comparison, if you call pmvnorm with the GenzBretz method
pmvnorm( ..., algorithm=GenzBretz)
the "error" attribute returns a sensible value, e.g. 3.855762e-6
You can get a positive answer with Miwa by cranking up the number of steps
it uses (to its maximum of 4098, for example)
pmvnorm( ..., algorithm=Miwa(steps=4098))
But this is still not a "correct" answer and the "error"
attribute is set
to NA to flag it.
HTH,
Eric
On Mon, Oct 2, 2017 at 6:42 PM, Eric Berger <ericjberger at gmail.com>
wrote:
> 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/posti
>> ng-guide.html
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
[[alternative HTML version deleted]]