[R] problem with a function

[li li]

Thanks very much for your reply.  The function is a bit long. I attached it.

rho.f () is in second text document. It uses rada1.mnorm() function which
is contained in the first document.

Thank you very very much!!!

2010/5/28 Dennis Murphy <djmuser at gmail.com>

> Hi:
>
> The problem is that input arguments such as corr[4] have to be evaluated
> within the body of your function, and apparently you haven't written it to
> do so. Unfortunately, I can't help further because my clairvoyance package
> is still in the concept development stage. In the meantime, it would be
> advisable if you could post a *minimal* version of the function that works
> with
> 0.3 but fails at corr[4] (as per instructions in the posting guide, which
> is linked
> at the bottom of this message. The simpler you make it for others to help,
> the more
> likely it will be that you'll get a satisfactory answer.
>
> HTH,
> Dennis
>
>  On Fri, May 28, 2010 at 6:52 AM, li li <hannah.hlx at gmail.com> wrote:
>
>>  Hi all,
>>   I have a function rho.f which gives a list of estimators. I have the
>> following problems.
>> rho.f(0.3) gives me the right answer. However, if I use rho.f(corr[4])
>> give
>> me a different
>> answer, even though corr[4]==0.3.
>>   This prevents me from using a for loop. Can someone give me some help?
>>    Thank you very much in advance.
>>                                                   Hannah
>>
>> > rho.f(0.3)
>> $est.1
>> [1] 0 0 0 0 0 0
>> $est.2
>> [1] 0 0 0 0 0 0
>> $est.3
>> [1] 0 0 0 0 0 0
>> $est.4
>> [1] 0 0 0 0 0 0
>> $est.5
>> [1] 0 0 0 0 0 0
>>
>> > corr <- seq(0,0.9, by=0.1)
>> > corr[4]
>> [1] 0.3
>>
>> > rho.f(corr[4])
>> $est.1
>> [1] 0.0 0.0 2.5 0.0 5.0 0.0
>> $est.2
>> [1] 0.00000 0.00000 0.00000 0.00000 3.72678 0.00000
>> $est.3
>> [1] 0.000000 0.000000 0.000000 0.000000 2.777778 0.000000
>> $est.4
>> [1]  0.00000  0.00000  0.00000  0.00000 13.88889  0.00000
>> $est.5
>> [1]  0.00000  0.00000  0.00000  0.00000 13.88889  0.00000
>> >
>>
>>        [[alternative HTML version deleted]]
>>
>>
>> ______________________________________________
>> 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<http://www.r-project.org/posting-guide.html>
>> and provide commented, minimal, self-contained, reproducible code.
>>
>
>
-------------- next part --------------
rdata1.mnorm<-function(m,n,mzero,mu0, mu1,rho )
{

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



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

J <- rep(1, m)

var.f <- function(rho) {
   (1-rho)*diag(m)+rho*J%*%t(J)
     }


set.seed(103)
x <- mvrnorm(n,mean, var.f(rho))

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




















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



###generate data

result <- rdata1.mnorm(m=10,n=6,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)}
        w4 <- sum(w)
         y <- list(vec=w,  sum=w4)  
         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)
}   

rho <- 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
                      }