[R] problem with a function

[li li]

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!!!


> rho.f(0.3)
$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

> corr <- seq(0,0.9, by=0.1)
> corr[4]
[1] 0.3
> rho.f(corr[4])
$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>

>
> On May 28, 2010, at 12:03 PM, Sarah Goslee wrote:
>
>  From your code:
>>>
>>
>> mean <- c(rep(mu0, mzero), rep(mu1,m-mzero))
>>
>> mean() is a function. If you overwrite it with data, you may mess
>> other things up -
>> any function you call that calls mean will now fail (at best).
>>
>
> Actually, it's bad but not quite that bad.
>
> > mean <- c(1,2,3,4)
> > mean(1:10)
> [1] 5.5
> > mean(mean)
> [1] 2.5
> > mean[3]
> [1] 3
> > mean
> [1] 1 2 3 4
>
> > apply(matrix(1:100, ncol=10),1, mean)
>  [1] 46 47 48 49 50 51 52 53 54 55
> > rm(mean)  # the interpreter "knows" to only remove the vector object
> > mean
> function (x, ...)
> UseMethod("mean")
> <environment: namespace:base>
> > rm(mean)  # and will refuse to remove the base function
> Warning message:
> In rm(mean) : object 'mean' not found
>
>
> Generally the interpreter can tell when a name is being intended as a
> function. Certainly when () follows the name, a function will be sought and
> other objects with identical names ignored.There are exceptions to that
> statement and your point is very well taken, but the main level of confusion
> is in the human brain rather than the R-interpreter. There used to be more
> partial name matching, but my reading of the NEWS items makes me think there
> is a shift away from that facility.  Other functions that people often
> mis-use as object names, generally without obvious deleterious effects:
>
> df  # the density of the F distribution
> c
> data
> sd
> var
> names
>
>
>> var.f <- function(rho) {
>>  (1-rho)*diag(m)+rho*J%*%t(J)
>> }
>>
>> var.f() is a complete function, except that m and J are not passed
>> as arguments. Instead, you rely on them being present in the
>> calling environment, and that is both dangerous and bad practice.
>>
>> Sarah
>>
>> On Fri, May 28, 2010 at 12:00 PM, li li <hannah.hlx at gmail.com> wrote:
>>
>>> I am not sure about "overwrite mean() with data".  My purpose was
>>> to generate random numbers that are from a multivariate normal
>>> distribution with the mean vector.
>>>
>>> For the var.f function, since I already specify m and J, so the only
>>> variable is really rho, so I wrote it as a function of rho only.
>>>
>>> Could you be a little more specific? Thanks a lot again.
>>> 2010/5/28 Sarah Goslee <sarah.goslee at gmail.com>
>>>
>>>>
>>>> There are a bunch of problems in your code:
>>>> you overwrite mean() with data, and that could screw things up.
>>>> you have a function var.f that isn't passed all the arguments it needs.
>>>> est.4 is defined several times, each overwriting the previous.
>>>>
>>>> First you need to clean up these sorts of problems since they can
>>>> lead to all kinds of bizarre results.
>>>>
>>>> Then, if you are still getting unexpected results, please send the list
>>>> a minimal example so that we can take a look.
>>>>
>>>> Sarah
>>>>
>>>
>> --
>> Sarah Goslee
>> http://www.functionaldiversity.org
>>
>> ______________________________________________
>> 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.
>>
>
> David Winsemius, MD
> West Hartford, CT
>
>
-------------- 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
                      }