From: Joris Meys <jorismeys_at_gmail.com>

Date: Fri, 18 Jun 2010 19:41:23 +0200

Date: Fri, 18 Jun 2010 19:41:23 +0200

See below.

On Fri, Jun 18, 2010 at 7:11 PM, li li <hannah.hlx_at_gmail.com> wrote:

> Dear all,

*> I am trying to calculate certain critical values from bivariate normal
**> distribution (please see the
**> function below).
**>
**> m <- 10
**> rho <- 0.1
**> k <- 2
**> alpha <- 0.05
**> ## calculate critical constants
**> cc_z <- numeric(m)
**> var <- matrix(c(1,rho,rho,1), nrow=2, ncol=2, byrow=T)
**> for (i in 1:m){
**> if (i <= k) {cc_z[i] <- qmvnorm((k*(k-1))/(m*(m-1))*alpha, tail="upper",
**> sigma=var)$quantile} else
**> {cc_z[i] <- qmvnorm((k*(k-1))/((m-i+k)*(m-i+k-1))*alpha,
**> tail="upper", sigma=var)$quantile}
**> }
**>
**>
**>
**> After the critical constants cc_z is calculated, I wanted to check whether
**> they are correct.
**>
**>
**>> ##check whether cc_z is correct
**>> pmvnorm(lower=c(cc_z[1], cc_z[1]),
**> upper=Inf,sigma=var)-(k*(k-1))/(n*(n-1))
*

Shouldn't this be

> pmvnorm(lower=c(cc_z[1], cc_z[1]),

+ upper=Inf,sigma=var)-(k*(k-1))/(m*(m-1))*alpha
[1] -5.87e-09

attr(,"error")

[1] 1e-15

attr(,"msg")

[1] "Normal Completion"

This still gives a bit of an error, but you have to take into account as well that the underlying algorithms use randomized quasi-MC methods, and that floating point issues can play here as well. So it looks to me that your calculations are correct.

Cheers

Joris

-- Joris Meys Statistical consultant Ghent University Faculty of Bioscience Engineering Department of Applied mathematics, biometrics and process control tel : +32 9 264 59 87 Joris.Meys_at_Ugent.be ------------------------------- Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php ______________________________________________ R-help_at_r-project.org mailing list 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.Received on Fri 18 Jun 2010 - 17:43:16 GMT

Archive maintained by Robert King, hosted by
the discipline of
statistics at the
University of Newcastle,
Australia.

Archive generated by hypermail 2.2.0, at Fri 18 Jun 2010 - 18:10:33 GMT.

*
Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-help.
Please read the posting
guide before posting to the list.
*