MCMChybridGP (version 5.4)

hybrid.sample: Sampling phase applying results from Exploratory phase.

Description

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.

Usage

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)

Arguments

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.

Value

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.

Details

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

References

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

See Also

hybrid.explore

Examples

# NOT RUN {
    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))
    }

# }
# NOT RUN {
    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)
# }