R help - Nov 2011 - Error in random walk Metroplis-hasting

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]]

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

Your code produces likelihood estimates of zero and divides them.
Hence NaN. You should decide for yourself whether this is a data
problem or an algorithmic problem.

In general try using options(error = recover) to debug when you get in
a situation like this: it's super helpful.

I'll happily give you more help if you make a good faith effort to do
what I suggested in my first reply.

Michael

On Wed, Nov 16, 2011 at 1:27 PM, Gyanendra Pokharel
<gyanendra.pokharel at gmail.com> wrote:
> Hi Michael thanks for the response. Following is my data set. Could you
> please see the code through this data set.
> ??? X indnum? x? y inftime
> 1???? 1????? 1? 1? 1?????? 0
> 2???? 2????? 2? 2? 1?????? 0
> 3???? 3????? 3? 3? 1?????? 0
> 4???? 4????? 4? 4? 1?????? 0
> 5???? 5????? 5? 5? 1?????? 0
> 6???? 6????? 6? 6? 1?????? 0
> 7???? 7????? 7? 7? 1?????? 0
> 8???? 8????? 8? 8? 1?????? 0
> 9???? 9????? 9? 9? 1?????? 0
> 10?? 10???? 10 10? 1?????? 0
> 11?? 11???? 11? 1? 2?????? 0
> 12?? 12???? 12? 2? 2?????? 0
> 13?? 13???? 13? 3? 2?????? 0
> 14?? 14???? 14? 4? 2?????? 0
> 15?? 15???? 15? 5? 2?????? 0
> 16?? 16???? 16? 6? 2?????? 0
> 17?? 17???? 17? 7? 2?????? 0
> 18?? 18???? 18? 8? 2?????? 0
> 19?? 19???? 19? 9? 2?????? 0
> 20?? 20???? 20 10? 2?????? 0
> 21?? 21???? 21? 1? 3?????? 0
> 22?? 22???? 22? 2? 3?????? 0
> 23?? 23???? 23? 3? 3?????? 0
> 24?? 24???? 24? 4? 3?????? 0
> 25?? 25???? 25? 5? 3?????? 0
> 26?? 26???? 26? 6? 3?????? 0
> 27?? 27???? 27? 7? 3?????? 0
> 28?? 28???? 28? 8? 3?????? 0
> 29?? 29???? 29? 9? 3?????? 0
> 30?? 30???? 30 10? 3?????? 0
> 31?? 31???? 31? 1? 4?????? 0
> 32?? 32???? 32? 2? 4?????? 0
> 33?? 33???? 33? 3? 4?????? 0
> 34?? 34???? 34? 4? 4?????? 0
> 35?? 35???? 35? 5? 4?????? 0
> 36?? 36???? 36? 6? 4?????? 0
> 37?? 37???? 37? 7? 4?????? 0
> 38?? 38???? 38? 8? 4?????? 0
> 39?? 39???? 39? 9? 4?????? 0
> 40?? 40???? 40 10? 4?????? 0
> 41?? 41???? 41? 1? 5?????? 0
> 42?? 42???? 42? 2? 5?????? 0
> 43?? 43???? 43? 3? 5?????? 2
> 44?? 44???? 44? 4? 5?????? 0
> 45?? 45???? 45? 5? 5?????? 0
> 46?? 46???? 46? 6? 5?????? 0
> 47?? 47???? 47? 7? 5?????? 0
> 48?? 48???? 48? 8? 5?????? 0
> 49?? 49???? 49? 9? 5?????? 0
> 50?? 50???? 50 10? 5?????? 0
> 51?? 51???? 51? 1? 6?????? 0
> 52?? 52???? 52? 2? 6?????? 0
> 53?? 53???? 53? 3? 6?????? 0
> 54?? 54???? 54? 4? 6?????? 0
> 55?? 55???? 55? 5? 6?????? 0
> 56?? 56???? 56? 6? 6?????? 0
> 57?? 57???? 57? 7? 6?????? 0
> 58?? 58???? 58? 8? 6?????? 0
> 59?? 59???? 59? 9? 6?????? 2
> 60?? 60???? 60 10? 6?????? 0
> 61?? 61???? 61? 1? 7?????? 0
> 62?? 62???? 62? 2? 7?????? 0
> 63?? 63???? 63? 3? 7?????? 0
> 64?? 64???? 64? 4? 7?????? 0
> 65?? 65???? 65? 5? 7?????? 0
> 66?? 66???? 66? 6? 7?????? 2
> 67?? 67???? 67? 7? 7?????? 0
> 68?? 68???? 68? 8? 7?????? 0
> 69?? 69???? 69? 9? 7?????? 0
> 70?? 70???? 70 10? 7?????? 0
> 71?? 71???? 71? 1? 8?????? 0
> 72?? 72???? 72? 2? 8?????? 2
> 73?? 73???? 73? 3? 8?????? 0
> 74?? 74???? 74? 4? 8?????? 2
> 75?? 75???? 75? 5? 8?????? 0
> 76?? 76???? 76? 6? 8?????? 2
> 77?? 77???? 77? 7? 8?????? 2
> 78?? 78???? 78? 8? 8?????? 0
> 79?? 79???? 79? 9? 8?????? 2
> 80?? 80???? 80 10? 8?????? 0
> 81?? 81???? 81? 1? 9?????? 0
> 82?? 82???? 82? 2? 9?????? 0
> 83?? 83???? 83? 3? 9?????? 0
> 84?? 84???? 84? 4? 9?????? 0
> 85?? 85???? 85? 5? 9?????? 0
> 86?? 86???? 86? 6? 9?????? 2
> 87?? 87???? 87? 7? 9?????? 2
> 88?? 88???? 88? 8? 9?????? 2
> 89?? 89???? 89? 9? 9?????? 2
> 90?? 90???? 90 10? 9?????? 0
> 91?? 91???? 91? 1 10?????? 0
> 92?? 92???? 92? 2 10?????? 0
> 93?? 93???? 93? 3 10?????? 0
> 94?? 94???? 94? 4 10?????? 0
> 95?? 95???? 95? 5 10?????? 2
> 96?? 96???? 96? 6 10?????? 2
> 97?? 97???? 97? 7 10?????? 1
> 98?? 98???? 98? 8 10?????? 2
> 99?? 99???? 99? 9 10?????? 2
> 100 100??? 100 10 10?????? 2
>
>
>
>
> Thanks
>
> On Wed, Nov 16, 2011 at 1:13 PM, R. Michael Weylandt
> <michael.weylandt at gmail.com> wrote:
>>
>> 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.
>> >
>
>