0

I am trying to write a function that can calculate Approximate Bayesian Computation using the Population Monte Carlo method. However, I ran into some troubles with my R code with the following error.

>head(ABC_PMC(2,5000,1,observedData,observedSummary))

Error in while (prior == 0) { : missing value where TRUE/FALSE needed
In addition: Warning message:
In log(importance.sample[2]) : NaNs produced

The code doesn't have a problem when I run it with one iteration, t=1. The issue lies with the second iteration and I think it is the assignment of importance.sammpling step. But I can't figure out what is the error as I tried to run the code manually and it works.

Any and all help will be appreciated.

# initialise variables:

    set.seed(12345)
    n=5000
    weight = rep(1,n)
    theta = matrix(NA, nrow = n, ncol = 2)
    for (i in 1:n) {
     fit = c(runif(1,-2,2), runif(1,0,4))
    theta[i,] = fit
    }
    
        
    #Calculate importance density for current t to use in next iteration
        old.prop.weight = weight / sum(weight)
        old.theta = theta
        Kernel.mean = cbind(old.theta[,1],old.theta[,2])
        Kernel.cov = cov(old.theta)   

Main loop:

      for (i in 1:n) {
            prior=0
            while (prior==0) {
              #Loop to test is prior=0

                choice = sample(1:n, size = 1, prob = old.prop.weight)
                importance.sample = rmvnorm(1,mean = Kernel.mean[choice,], sigma = Kernel.cov)
                #logsigma to ensure theta2 are all nonnegative
                if (importance.sample[2] != 0)
                  importance.sample[2] = log(importance.sample[2])
                fit = c(importance.sample[1],importance.sample[2])
              }
              
              prior = as.numeric(dunif(fit[1],min = -2,max = 2) * dunif(exp(fit[2]),min = 0,max = 4))
            }
            
           
          
          theta[i,] = fit
          ##End of step 2 loop
        }

Edit: Tried to make the code into a reproducible example.

janw
  • 8,758
  • 11
  • 40
  • 62
Stongals
  • 1
  • 2
  • (1) can we have a [reproducible example](http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example) please? (2) have you tried using `options(error=recover)` to drop you into the browser so you can see what's going on at the point where the error occurs? – Ben Bolker Nov 01 '15 at 15:14
  • @BenBolker I have tried using options(error=recover) and it brought me to the prior step. The error caused there is caused by the previous importance.sampling assignment step though – Stongals Nov 01 '15 at 15:32
  • thanks for trying to make it reproducible, but: (1) you need to define `n` and `theta`; (2) `rmvnorm` is from the `mvtnorm` package (you could also use `MASS::mvrnorm()`; (3) as written, the code will get stuck indefinitely in the `while (prior==0) { }` clause ... – Ben Bolker Nov 01 '15 at 15:38
  • Apologies. @BenBolker . n=5000, theta is defined in initialisation step to be a 5000X2 matrix – Stongals Nov 01 '15 at 16:13

0 Answers0