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

[Mike Lawrence]

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