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.
