Skip to content
Prev 222440 / 398500 Next

problem with a function

I modified my codes. However it looks like it still has the same problem.
Again, rho.f(0.3) gives the right answer. rho.f(corr[4])
gives wrong answer even though corr[4]==0.3.


The codes are attached.

Thank you very much!!!
$est.1
 [1] 0.0000000 0.0000000 0.0000000 0.0000000 0.8333333 0.0000000 0.0000000
 [8] 0.0000000 1.6666667 0.0000000
$est.2
 [1] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
 [9] 1.198658 0.000000
$est.3
 [1] 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000 0.000000
 [9] 0.862069 0.000000
$est.4
 [1]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
 [9] 12.93103  0.00000
$est.5
 [1]  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  0.00000
 [9] 12.93103  0.00000
[1] 0.3
$est.1
 [1] 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000 0.0000000
 [8] 0.0000000 0.8333333 0.0000000
$est.2
 [1] 0 0 0 0 0 0 0 0 0 0
$est.3
 [1] 0 0 0 0 0 0 0 0 0 0
$est.4
 [1] 0 0 0 0 0 0 0 0 0 0
$est.5
 [1] 0 0 0 0 0 0 0 0 0 0
2010/5/28 David Winsemius <dwinsemius at comcast.net>
-------------- next part --------------
rdata1.mnorm<-function(m,n,mzero,mu0, mu1,rho )
{

## ARGUEMENTS
 # n: sample size
 # m: dimension of multivariate normal



library(MASS)
 
mean1 <- c(rep(mu0, mzero), rep(mu1,m-mzero))


var.f <- function(rho, x,y) {
   (1-rho)*diag(x)+rho*y%*%t(y)
     }

J <- rep(1, m)

set.seed(103)
x <- mvrnorm(n,mean1, var.f(rho,x=m, y=J))

theta <- matrix(0, nrow=n, ncol=m)
for (i in 1:m){theta[,i]<- mean1[i]>1}
data<-list(s=theta, o=x)
return(data) 
                }




















-------------- next part --------------
 
rho.f <- function(rho){



###generate data

result <- rdata1.mnorm(m=30,n=10,mzero=5,mu0=0,mu1=2,rho=rho)
o <- result$o
s <- result$s

### the p-values
pv<-1-pnorm(o, 0, 1)

m <- dim(pv)[2]
n <- dim(pv)[1]


lambda <- 0.96



w1 <- apply(pv>lambda, 1, sum)

est.1 <-  w1/((1-lambda)*m)



w2 <- numeric(n)

for (i in 1:n){w2[i]<- choose(w1[i],2)}

est2 <- 2*w2/((1-lambda)^2*m*(m-1))

est.2 <- sqrt(est2)




est.3 <- numeric(n)

for (i in 1:n){
                if (est.1[i]>0)
                {est.3[i]<- est2[i]/est.1[i]}
                 else
                {est.3[i] <- 0}
                     }


est.4 <- numeric(n)



w4.f <- function(x, lambda){
        k <- length(x)-1
      w <- numeric(k)
      for (i in 1:k){
         w[i] <- (x[i]>lambda)*(x[i+1]>lambda)}
        s1 <- sum(w)
         y <- list(vec=w,  sum=s1)  
         return(y) 
             }

w4 <- numeric(n)
for (i in 1:n){ w4[i] <- as.numeric(w4.f(pv[i,], lambda)$sum)}

est4 <- w4/((1-lambda)^2*(m-1))

est.4 <- numeric(n)

for (i in 1:n){if (est.1[i]>0)
                {est.4[i]<- est4[i]/est.1[i]}
                 else est.4[i] <- 0}





st.pv <- t(apply(pv,1,sort))



w5 <- numeric(n)
for (i in 1:n){ w5[i] <- as.numeric(w4.f(st.pv[i,], lambda)$sum)}

est5 <- w5/((1-lambda)^2*(m-1))

est.5 <- numeric(n)

for (i in 1:n){if (est.1[i]>0)
                {est.5[i]<- est5[i]/est.1[i]}
                 else est.5[i] <- 0}


y <- list(est.1=est.1, est.2=est.2, est.3=est.3, est.4=est.4,
est.5=est.5)
return(y)
}   

corr <- seq(0,0.9, by=0.1)
 k<- length(rho)
 

est.1 <- matrix(0, nrow=k, ncol=n)


est.2 <- matrix(0, nrow=k, ncol=n)


est.3 <- matrix(0, nrow=k, ncol=n)


est.4 <- matrix(0, nrow=k, ncol=n)


est.5 <- matrix(0, nrow=k, ncol=n)

for (i in 1:k){
      result <- rho.f(rho[i])
                  est.1[i,] <- result$est.1
                  est.2[i,] <- result$est.2
                  est.3[i,] <- result$est.3
                  est.4[i,] <- result$est.4
                  est.5[i,] <- result$est.5
                      }