[R] help a newbie with a loop
Simon Knapp
sleepingwell at gmail.com
Mon Jul 3 15:10:55 CEST 2006
On 7/3/06, Boks, M.P.M. <M.P.M.Boks at umcutrecht.nl> wrote:
>
>
> Hi,
>
> I am new in R and stumbled on a problem my (more experienced) friends
> can not help with with. Why isnt this code working?
>
> The function is working, also with the loop and the graph appears,
>
> only when I build another loop around it (for different values of p) ,
> R stays in a loop?
>
> Can't it take more then 2 loops in one program?
You can have as many as you like (as far as I know)
>
>
> powerb<-function(x,sp2,a,b,b1,m)
> { sx<-(sum(x^2)-(sum(x)^2)/length(x))/length(x)
> n0<-ceiling((((qnorm(1-(a/2))+qnorm(1-b))/b1)^2)*(((m+1)/m)*sp2/sx))
> repeat
> {
> n1<-ceiling((((qt(1-(a/2),n0-4)+qt(1-b,n0-4))/b1)^2)*(((m+1)/m)*sp2/sx))
> if(n0==n1) break
> n0<-n1
> }
> return(c(sx,n1))
> }
>
> x<-rnorm(1000,0,1)
> x<-x[order(x)]
>
> res<-matrix(0,1000,2)
>
>
> #use the function and plot for different values of ind and p
> for ( p in c(0.05,0.10,0.15,0.20,0.25,0.30,0.40,0.50))
> { risk<-p*(2-p)
> nonrisk<-(1-p)^2
> m<-nonrisk/risk
>
> for (ind in 1:500)
> {res[ind,]<-powerb(x[c(1:(500-ind),(500+ind):1000)],4,0.05,0.20,0.1,m)}
>
> plot(res[,1],res[,2],type="p",ylab="n",xlab="var(x)",main="b=0.1,power=0
> .80,alpha=0.05,dominant met p=0.25")}
>
I modified your function as follows:
powerb<-function(x,sp2,a,b,b1,m){
sx<-(sum(x^2)-(sum(x)^2)/length(x))/length(x)
n0<-ceiling((((qnorm(1-(a/2))+qnorm(1-b))/b1)^2)*(((m+1)/m)*sp2/sx))
repeat{
n1<-ceiling((((qt(1-(a/2),n0-4)+qt(1-b,n0-4))/b1)^2)*(((m+1)/m)*sp2/sx))
if(n0==n1) break else cat("p = ", p, ", ind = ", ind, ", n0 =
", n0, ", n1 = ", n1, "\n", sep='')
n0<-n1
}
return(c(sx,n1))
}
and got, for example, the following output:
....
p = 0.2, ind = 278, n0 = 2288, n1 = 2289
p = 0.2, ind = 278, n0 = 2289, n1 = 2288
p = 0.2, ind = 278, n0 = 2288, n1 = 2289
p = 0.2, ind = 278, n0 = 2289, n1 = 2288
p = 0.2, ind = 278, n0 = 2288, n1 = 2289
p = 0.2, ind = 278, n0 = 2289, n1 = 2288
p = 0.2, ind = 278, n0 = 2288, n1 = 2289
p = 0.2, ind = 278, n0 = 2289, n1 = 2288
...
so your in an infinite loop. Maybe you can check the difference
between n0 and n1 instead, or perhaps better, do something like:
powerb<-function(x,sp2,a,b,b1,m){
sx<-(sum(x^2)-(sum(x)^2)/length(x))/length(x)
n0<-ceiling((((qnorm(1-(a/2))+qnorm(1-b))/b1)^2)*(((m+1)/m)*sp2/sx))
n0.1 <- n1.1 <- Inf
repeat{
n1<-ceiling((((qt(1-(a/2),n0-4)+qt(1-b,n0-4))/b1)^2)*(((m+1)/m)*sp2/sx))
if(n0==n1) break
if(n0.1==n1 && n1.1==n0) break
n0.1 <- n0
n1.1 <- n1
n0<-n1
}
return(c(sx,n1))
}
More information about the R-help
mailing list