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))
[1] 0.001111025
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"> pmvnorm(lower=c(cc_z[3], cc_z[3]),
upper=Inf,sigma=var)-(k*(k-1))/((n-1)*(n-2))
[1] 0.001388974
attr(,"error")
[1] 1e-15
attr(,"msg")
[1] "Normal Completion"
As you can see, there is some error. Supposedly, the answer should be zero
for all cc_z[i]
in the above checking function.
Is there something wrong with my code? Can someone give some help?
Thank you very much!
Hannah
[[alternative HTML version deleted]]
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 : helpdesk.ugent.be/e-maildisclaimer.php