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 <statsstud...@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@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. ~
______________________________________________
R-help@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.
