[R] Error in random walk Metroplis-hasting
R. Michael Weylandt
michael.weylandt at gmail.com
Thu Nov 17 07:18:50 CET 2011
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.
>> >
>
>
More information about the R-help
mailing list