R help - Apr 2009 - "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

_________________________________________________________________


	[[alternative HTML version deleted]]

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

Your problem is that with the alpha & beta you've specified

(((alpha-1)^(alpha-1))*((beta-1)^(beta-1)))/((beta+alpha-2)^(alpha+beta-2))

is

Inf/Inf

which is NaN.



On Sun, Apr 12, 2009 at 5:39 PM, Mary Winter <statsstudent at hotmail.com>
wrote:
>
>
> ?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
>
> _________________________________________________________________
>
>
> ? ? ? ?[[alternative HTML version deleted]]
>
> ______________________________________________
> R-help at 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.
>


-- 
Mike Lawrence
Graduate Student
Department of Psychology
Dalhousie University

Looking to arrange a meeting? Check my public calendar:
http://tinyurl.com/mikes-public-calendar

~ Certainty is folly... I think. ~