[R] slice sampler

[Tong Wang]

Hi there,  
     I wrote a slice sampler (one dimension) ,  After testing it, I found it working well on some standard distributions, but problematic with others.  I have thought about it but really can't figure out how to improve it, so I attached the code here
and hope that someone here is willing to try it out and give me some feedback. I would really appreciate it. 
     I basicly followed neal 03 to write the code. The code might looks a little strange as I allowed sampling a vector of parameters  ,   that is , it could update multiple targets with different parameter but same form of distribution. A example
of usage is as follows:    
      output[i] <- slice1d(dnorm(x,1,1,log=TRUE), x0=output[i-1],w=2,m=6,trunc=c(-Inf,Inf),prec=1e-4) 
 things to notice here are: expression has to be on LOG scale, w is the expanding wideth, m is the # of expanding trials,
prec is the precision, this is used to control against those distributions with huge kurtosis (with very sharp peak),  so that the sampler wouldn't spend a lot of time approaching the right value. 

One example that it doesn't work well on is    dgamma(x,.1,.1,log=TRUE)

Thanks a lot in advance. 
best,

CODE: 

slice1d <- function (expr,x0=NULL,w,m,trunc=c(-Inf,Inf),prec=.0001){
    sexpr <- substitute(expr)
    if (!(is.call(sexpr) && match("x", all.vars(sexpr), nomatch = 0))) 
      stop("'expr' must be a function or an expression containing 'x'")
    expr <- sexpr
    
    y.log <- eval(expr, envir=list(x=x0),enclos = if(is.null(list(...)$env)) parent.frame() else list(...)$env)-rexp(1)
    k <- length(y.log)
    if(k!=length(x0))
      stop("The dimension of y.log and x0 have to agree")

    U <- runif(k)
    L <- x0-w*U
    R <- L+w
    U <- runif(k)
    J <- floor(m*U)
    K <- m-1-J

    temp <- (L<trunc[1])   
    L <- as.vector(L)
    L[temp] <- trunc[1]
    vec.cond <- rep(!temp,k)
    
    while(sum(vec.cond)!=0)
      {
        vec.cond <- ((J>0)&(y.log<eval(expr,list(x=L),enclos = if(is.null(list(...)$env)) parent.frame() else list(...)$env)))
        L[vec.cond] <- L[vec.cond]-w
        temp <- (L<trunc[1])
        L[temp] <- trunc[1]
        J[temp] <- 0
        J[vec.cond] <- J[vec.cond]-1}

    temp <- (R>trunc[2])   
    R <- as.vector(R)
    R[temp] <- trunc[2]
    vec.cond <- rep(!temp,k)

    while(sum(vec.cond)!=0)
      {vec.cond <- ((K>0)&(y.log<eval(expr,list(x=R),enclos = if(is.null(list(...)$env)) parent.frame() else list(...)$env)))
       R[vec.cond] <- R[vec.cond]+w
       temp <- (R>trunc[2])
       R[temp] <- trunc[2]
       K[temp] <- 0
       K[vec.cond] <- K[vec.cond]-1}

    vec.cond <- rep(FALSE,k)
    x1 <- rep(NA,k)
    count <- 0
    repeat{
      U <- runif((k-sum(vec.cond)))
      x1[!vec.cond] <- L[!vec.cond]+U*(R[!vec.cond]-L[!vec.cond])
      vec.cond <- ((y.log<eval(expr,list(x=x1),enclos = if(is.null(list(...)$env)) parent.frame() else list(...)$env))|(abs(x1-x0)<prec)) 
      if(sum(vec.cond)==k) {return(x1)}
      L[x1<=x0] <- x1[x1<=x0]
      R[x1>x0] <- x1[x1>x0]
      count <- count+1
      if (count>=10000) stop("dead loop")
    }
  }