"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



Get the New Internet Explore 8 Optimised for MSN. Download Now