[R] Error in random walk Metroplis-hasting

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