-----Original Message-----
From: r-help-bounces at r-project.org
[mailto:r-help-bounces at r-project.org] On Behalf Of Mary Winter
Sent: Sunday, April 12, 2009 1:39 PM
To: r-help at 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