An embedded and charset-unspecified text was scrubbed... Name: not available URL: <https://stat.ethz.ch/pipermail/r-help/attachments/20111116/fb859f12/attachment.pl>
Error in random walk Metroplis-hasting
2 messages · Gyanendra Pokharel, R. Michael Weylandt
Three things: 1) This isn't reproducible without your data file. Work out a minimal reproducible example -- I bet you'll find your answer along the way -- and supply it. 2) What do the warnings say: they are usually pretty good at directing you to trouble. 3) Don't use attach. Seriously -- just load the data in once and pass it around where needed. (It's going to slow you down to re-load it each time (~400 times) you call likelihood.) Michael PS -- As a style point, might I suggest you put more spaces around the assignment operator "<-" it's hard (for both a human and the R interpreter) to tell is x<-3 means to set x equal to 3 or to test if x is less than "-3" and sometimes this leads to very tricky bugs. On Wed, Nov 16, 2011 at 10:36 AM, Gyanendra Pokharel
<gyanendra.pokharel at gmail.com> wrote:
Hi R community,
I have some data set and construct the likelihood as follows
likelihood <- function(alpha,beta){
? ?lh<-1
? ?d<-0
? ?p<-0
? ?k<-NULL
? ?data<-read.table("epidemic.txt",header = TRUE)
? ?attach(data, warn.conflicts = F)
? ?k <-which(inftime==1)
? ?d <- (sqrt((x-x[k])^2+(y-y[k])^2))^(-beta)
? ?p<-1 - exp(-alpha*d)
? ?for(i in 1:100){
? ? ? ?if(i!=k){
? ? ? ? ? ?if(inftime[i]==0){
? ? ? ? ? ? ? ?lh<-lh*(1-p[i])
? ? ? ? ? ?}
? ? ? ? ? ?if(inftime[i]==2){
? ? ? ? ? ? ? ?lh<-lh*p[i]
? ? ? ? ? ?}
? ? ? ?}
? ?}
? ?return(lh)
}
?Then, I want to simulate this by using random walk Metropolis- Hasting
algorithm with the single parameter update. i have the following R function
mh.epidemic <-function(m,alpha, beta){
? ? ?z<- array(0,c(m+1, 2))
? ?z[1,1] <- alpha
? ?z[1,2] <- beta
? ?for(t in 2:m){
? ? ? ?u <- runif(1)
? ? ? ?y1 <- rexp(1,z[t-1,1])
? ? ? ?y2 <-rexp(1,z[t-1,2])
? ? ? ?z[t,1] <- likelihood(y1,z[t-1,2])/likelihood(z[t-1,1],z[t-1,2])
? ? ? ?a1 <-min(1,z)
? ? ? ?if(u<=a1) z[t,] <- y1 else
? ? ? ?z[t,2] <-z[t-1,1]
? ? ? ?z[t,2]<-likelihood(z[t,1], y2)/likelihood(z[t,1],z[t-1,2])
? ? ? ?accept <-min(1,z)
? ? ? ?if(u<=accept) z[t,] <- y2 else
? ? ? ?z[t,2] <-z[t-1,1]
? ?}
? ?z
}
mh.epidemic(100, 1,2)
when I run it shows the following error:
Error in if (u <= accept) z[t, ] <- y2 else z[t, 2] <- z[t - 1, 1] :
?missing value where TRUE/FALSE needed
I know it is due to the NaN produce some where. So I tried by running the
code separately, as follows
m<- 100
? ? alpha <-1
? ? beta <- 2
? ? z<- array(0,c(m+1, 2))
? ? z[1,1] <- alpha
? ? z[1,2] <- beta
? ? for(t in 2:m){
+ ? ? ? ? u <- runif(1) + ? ? ? ? y1 <- rexp(1,z[t-1,1]) + ? ? ? ? y2 <-rexp(1,z[t-1,2]) + ? ? ? ? z[t,1] <- likelihood(y1,z[t-1,2])/likelihood(z[t-1,1],z[t-1,2]) + ? ? ? ? a1 <-min(1,z) + } There were 50 or more warnings (use warnings() to see the first 50)
y1
[1] NaN
y2
[1] NaN why, this simulation from exponential function is produce NaN? If some one help me it will be great. ? ? ? ?[[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.