hybrid.sample: Sampling phase applying results from Exploratory phase.
In MCMChybridGP: Hybrid Markov Chain Monte Carlo using Gaussian Processes

Description Usage Arguments Details Value Note Author(s) References See Also Examples

Sampling phase in hybrid MCMC, which takes the output from hybrid.explore to samples from the target distribution supplied. The number of chains, Leapfrog moves and Gaussian process parameters are the same as used in hybrid.explore, or the values may be updated here. For a target distribution time consuming to evaluate, the target can be replaced completely by the Gaussian process approximation in some or all of the chains. Bounds supplied act as reflecting barriers.

hybrid.sample(Explore, n = 1000, replace.target = c(0, 1, 2),
              lb = NULL, ub = NULL, L = NULL, delta = NULL,
              nchains = NULL, T.mult = NULL, maxleap = NULL,
              r = 5, nswaps = NULL, graph = FALSE)

`Explore`	Output from `hybrid.explore`. A list object consisting of `X`, `y`, `f` and `GPfit`. `X` is the set of points from which a Gaussian process fit is gained. Then `y` is the corresponding values of the target log-density function `f`, with further arguments `...` as supplied to `hybrid.explore`. If `GPfit` is `NULL` then a new Gaussian process fit is determined.
`n`	The number of sampling iterations.
`L`	An optional integer argument passed from `hybrid.explore`. The number of steps used in Leapfrog moves.
`delta`	An optional numerical argument passed from `hybrid.explore`. The size of steps used in Leapfrog moves.
`lb`	An optional numeric argument passed from `hybrid.explore`. Lower (finite) bounds placed on X. (If supplied, then `ub` must also be supplied and finite.)
`ub`	An optional numeric argument passed from `hybrid.explore`. Upper (finite) bounds placed on X. (If supplied, then `lb` must also be supplied and finite.)
`nchains`	An optional integer argument passed from `hybrid.explore`. The number of MCMC parallel chains to be used.
`T.mult`	An optional integer argument passed from `hybrid.explore`. The number of Parallel Tempering chains to be used.
`maxleap`	An optional numerical argument passed from `hybrid.explore`. Gives the maximum standard deviation of `f(x)` that a point can vary from points in Gaussian process fit.
`r`	An optional numerical argument (default 5). A penalty factor on points straying from the region of Gaussian process fit, when the target distribution is replaced.
`nswaps`	An optional integer argument passed from `hybrid.explore`. The number of repititions of swaps proposed between MCMC chains.
`replace.target`	The sampling scheme to be used (0, 1 or 2) in acceptance of MCMC proposals. Where 0 represents using the true target distribution in all chains. 1 (default) represents using the true target distribution only in the primary chain (having temperature 1). 2 represents replacing the target distribution in all chains by the Gaussian process approximation.
`graph`	An optional boolean argument (default is FALSE). Request graphical progress display during the sample phase.

The method used in hybrid.sample is described in Fielding, Nott and Liong (2011).

A list is returned consisting of the following.

`SAMP`	A matrix with rows corresponding to sampled points generated from the target distribution.
`y`	A column of the corrresponding values of the log-density of the target distribution.
`acceptance`	A column of 0 (rejected) and 1 (accepted) giving a record of sampling proposal acceptance.
`function.calls`	The number of function calls to evaluate the true log-density.

A record is kept throughout a run of hybrid.sample stored as a global variable list, hybrid.sample.out. Useful for a run stopped prematurely.

The method used in hybrid.sample gives extensions to the work of Rasmussen (2003) and is described in Fielding, Nott and Liong (2011).

Mark J. Fielding <mark.fielding@gmx.com>

"Efficient MCMC Schemes for Computationally Expensive Posterior Distributions", Fielding, Nott and Liong (2011).

hybrid.explore

    mu1 <- c(-1, -1)
    mu2 <- c(+1, +1)
    sigma.sq <- 0.16
    ub <- c(1.5, 3)
    X0 <- matrix(c(-2,-1, 0,-2, 0, 1, 0, 1, 1,
                   -2,-1,-2, 0, 0, 0, 2, 1, 2), ncol = 2)
    f <- function(x) {
        px <- 1/4/pi/sqrt(sigma.sq) * exp(-1/2/sigma.sq *
            sum((x - mu1)^2)) + 1/4/pi/sqrt(sigma.sq) *
            exp(-1/2/sigma.sq * sum((x - mu2)^2))
        return(log(px))
    }


    explore.out <- hybrid.explore(f, X0, ub=ub, n=150, graph=TRUE)
    sample.out <- hybrid.sample(explore.out, n=500, graph=TRUE)

    opar <- par(mfrow=c(2,1))
    plot(density(sample.out$SAMP[,1]), xlab="x1", ylab="f(x)")
    plot(density(sample.out$SAMP[,2]), xlab="x2", ylab="f(x)")
    par(opar)