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 <dwinsem...@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....@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.gos...@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@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
>
>
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)
}
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
}
______________________________________________
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.
