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
_________________________________________________________________
[[alternative HTML version deleted]]
