> -----Original Message----- > From: r-help-boun...@r-project.org > [mailto:r-help-boun...@r-project.org] On Behalf Of Mary Winter > Sent: Sunday, April 12, 2009 1:39 PM > To: r-help@r-project.org > Subject: [R] "taking the log then later exponentiate the result" query > > > > Hi, > > I am trying to figure out the observed acceptance rate and M, > using generalised rejection sampling to generate a sample > from the posterior distribution for p. > > I have been told my code doesn't work because I need to > "take the log of the expression for M, evaluate it and then > exponentiate the result." This is because R is unable to > calculate high powers such as 545.501. > > As you can see in my code I have tried to taking the log of M > and then the exponential of the result, but I clearly must be > doing something wrong. > I keep getting the error message: > > Error in if (U <= ratio/exp(M)) { : missing value where > TRUE/FALSE needed > > Any ideas how I go about correctly taking the log and then > the exponential? > > rvonmises.norm <- function(n,alpha,beta) { > out <- rep(0,n) > counter <- 0 > total.sim <- 0 > p<-alpha/(alpha+beta) > M > <-log((((alpha-1)^(alpha-1))*((beta-1)^(beta-1)))/((beta+alpha > -2)^(alpha+beta-2))) > while(counter < n) { > total.sim <- total.sim+1 > proposal <- runif(1) > if(proposal >= 0 & proposal <= 1) { > U <- runif(1) > ratio<- (p^(alpha-1))*((1-p)^(beta-1)) > if(U <=ratio/exp(M)) { > counter <- counter+1 > out[counter] <- proposal > } > } > } > obs.acc.rate <- n/total.sim > return(out,obs.acc.rate,M) > } > set.seed(220) > temp <- rvonmises.norm(10000,439.544,545.501) > print(temp$obs.acc.rate) > > Louisa > >
I think when someone told you to take the log of the calculation, they meant for you to simplify the logarithmic calculation algebraically so that you are not exponentiating large numbers. Try changing your calculation of M to (I think this right) M <- (alpha-1)*log(alpha-1) + (beta-1)*log(beta-1) - (alpha+beta-2)*log(alpha+beta-2) Hope this is helpful, Dan Daniel Nordlund Bothell, WA USA ______________________________________________ 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.