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