[R] "taking the log then later exponentiate the result" query

[Daniel Nordlund]

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

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