[R] (Fwd) Re: priority of operators in the FOR ( ) statement

[Petr Pikal]

Hi 

On 23 Aug 2005 at 12:03, Ravi.Vishnu at outokumpu.com wrote: 

> Dear All, 
> I spent an entire evening in debugging a small, fairly simple 
program 
> in R - without success. It was my Guru in Bayesian Analysis, 
Thomas 
> Fridtjof, who was able to diagonose the problem. He said that it 
took 
> a long time for him also to locate the problem. This program 
> illustrates in some ways the shortcomings of the error messages 
that R 

<snip>

> 
*****************************************************
***************** 
> ***' # Bayesian Data Analysis ## 
> source("M:/programming/Rfolder/Assignments/fortest.txt") 
>  
> # #Remove all objects from the workspace 
> rm(list=ls()) 
> # #We will also try to note the time that the program takes 
> # #We will start the clock at starttime 
> starttime <- proc.time()[3]; 
>  

> my.function<-function(x) { 
> s2<-sqrt(2); 
> if ((x>=0) & (x<s2)) return(x/2) 
> else  
> if ((x>=s2) & (x<1+s2)) return(0.2) 
> else  
> if ((x>=1+s2) & (x<1.5+s2)) return(0.6) 
> else  
> if ((x>1.5+s2) | (x<0)) return(0) 
> } 

I also wonder if this function computes what is intended 

maybe this 

vec<-c(0,sqrt(2),sqrt(2)+1,sqrt(2)+1.5) 

myfun<-function(x,vec) { 

y<-x/2 
y0<-c(0,1,.2,.6,0)[findInterval(x,vec)+1] 
pos<-which(y0%in%1) 
y0[pos]<-y[pos] 
y0 
} 

vec<-c(0,sqrt(2),sqrt(2)+1,sqrt(2)+1.5) 

> x<-rnorm(1000000) 
> system.time(my1<-myfun(x,vec)) 
[1] 1.62 0.05 1.67   NA   NA 
> 

will do it more efficiently. 

HTH 
Petr 

>  
> alphayx<-function(y,x) { 
> fy<-my.function(y) 
> fx<-my.function(x) 
> fyx<-fy/fx 
> # to account for 0/0 division 
> if (is.na(fyx)) fyx<-0 
> #fyx<-ifelse(is.na(fyx),0,fyx); 
> alpha<-min(1,fyx) 
> return(alpha) 
> } 
>  
> sigma<-0.5; 
> #nr is the number of iterations 
> nr<-20 
> x<-numeric(nr); 
> x[1]<-1; 
> t<-1:nr; 
>  
> for (i in 1:nr-1) { 
> xi<-x[i]; 
> yi<-rnorm(1,mean=xi,sd=sigma); 
> ui<-runif(1,0,1); 
> ualphai<-alphayx(yi,xi); 
> xn<-ifelse(ui<=ualphai,yi,xi); 
> x[i+1]<-xn; 
> } 
>  
> plot(t,x,type="p") 
>  
> endtime<-proc.time()[3]; 
> elapsedTime<-endtime-starttime; 
> cat("Elapsed time is", elapsedTime, "seconds", "\n") 
> 
*****************************************************
> ______________________________________________ 
> R-help at stat.math.ethz.ch mailing list 
> https://stat.ethz.ch/mailman/listinfo/r-help 
> PLEASE do read the posting guide! 
> http://www.R-project.org/posting-guide.html 

Petr Pikal
petr.pikal at precheza.cz