"taking the log then later exponentiate the result" query

-----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

?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