# What follows is the precise Splus code used for the Beta example in the
# research paper "Perfect Forward Simulation via Simulated Tempering",
# by S.P. Brooks, Y. Fan, and J.S. Rosenthal.


function(alpha = 25, beta = 75, n = 1)
{
        y <- array(NA, n)
        end <- 0
        mode <- (alpha-1)/(alpha+beta-2)
        c <- 1/dbeta(mode,alpha,beta)
        w1 <- 1/(c + 1)
        w0 <- 1 - w1
        astar <- min(1, (c * w1)/w0)
        epsilon <- astar/2
        for(j in 1:n) {
                count <- 0
                while(end != 1) {
                        count <- count + 1
                        time <- 1 + rgeom(1, epsilon)
                        tau <- 1
                        x <- runif(1)
                        i <-  - time + 1
                        while(i < 0) {
                                # Do tau update
                                acc <- 1
                                v0 <- 0
                                u <- runif(1)
                                if(w1 == 1/(c + 1)) {
                                        u <- 0
                                }
                                if(u <= 0.5) {
                                        v <- 0
                                }
                                else {
                                        v <- 1
                                }
                                if(v == 1) {
                                        if(astar < 1) {
                                                v0 <- runif(1, astar, 1)
                                        }
                                        else {
                                                v0 <- 1
                                        }
                                        if(tau == 0) {
                                                acc <- w1/(dbeta(x, alpha,
                                                        beta) * w0)
                                        }
                                }
                                else {
                                        v0 <- runif(1, 0, 1)
                                        if(tau == 1) {
                                                acc <- (dbeta(x, alpha, beta) *
                                                        w0)/w1
                                        }
                                }
                                acc <- min(1, acc)
                                if(v0 <= acc) {
                                        tau <- v
                                }
                                # Next update state
                                xnew <- runif(1)
                                if(tau == 0) {
                                        acc <- dbeta(xnew, alpha, beta)/dbeta(
                                                x, alpha, beta)
                                }
                                else {
                                        acc <- 1
                                }
                                if(runif(1) <= acc) {
                                        x <- xnew
                                }
                                i <- i + 1
                        }
                        if(tau == 0) {
                                end <- 1
                        }
                }
                y[j] <- x
                end <- 0
        }
        y
}