See below.
On Fri, Jun 18, 2010 at 7:11 PM, li li <hannah....@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.m...@ugent.be
-------------------------------
Disclaimer : http://helpdesk.ugent.be/e-maildisclaimer.php
______________________________________________
R-help@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.
