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

From: Petr Pikal <petr.pikal_at_precheza.cz>
Date: Wed 24 Aug 2005 - 21:25:36 EST


Hi

On 23 Aug 2005 at 12:03, Ravi.Vishnu@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@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@precheza.cz



R-help@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 Received on Wed Aug 24 21:29:19 2005

This archive was generated by hypermail 2.1.8 : Sun 23 Oct 2005 - 15:49:44 EST