makeCPMSampler {cPseudoMaRg}

correlated pseudo-marginal: generates functions that output a big vector

Description

correlated pseudo-marginal: generates functions that output a big vector

Usage

makeCPMSampler(
  paramKernSamp,
  logParamKernEval,
  logPriorEval,
  logLikeApproxEval,
  yData,
  numU,
  numIters,
  rho = 0.99,
  storeEvery = 1,
  nansInLLFatal = TRUE
)

Arguments

paramKernSamp

function(theta) -> theta proposal

logParamKernEval

function(oldTheta, newTheta) -> logDensity.

logPriorEval

function(theta) -> logDensity.

logLikeApproxEval

function(y, thetaProposal, uProposal) -> logApproxDensity.

yData

the observed data

numU

integer number of u samples

numIters

integer number of MCMC iterations

rho

correlation tuning parameter (-1,1)

storeEvery

increase this integer if you want to use thinning

nansInLLFatal

terminate the entire chain on NaNs, or simply disregard sample

Value

vector of theta samples

Examples

# sim data
realTheta1 <- .2 + .3
realTheta2 <- .2
realParams <- c(realTheta1, realTheta2)
numObs <- 10
realX <- rnorm(numObs, mean = 0, sd = sqrt(realTheta2))
realY <- rnorm(numObs, mean = realX, sd = sqrt(realTheta1 - realTheta2))
# tuning params
numImportanceSamps <- 1000
numMCMCIters <- 1000
randomWalkScale <- 1.5
recordEveryTh <- 1
sampler <- makeCPMSampler(
 paramKernSamp = function(params){
   return(params + rnorm(2)*randomWalkScale)
 },
 logParamKernEval = function(oldTheta, newTheta){
   dnorm(newTheta[1], oldTheta[1], sd = randomWalkScale, log = TRUE)
   + dnorm(newTheta[2], oldTheta[2], sd = randomWalkScale, log = TRUE)
 },
 logPriorEval = function(theta){
   if( (theta[1] > theta[2]) & all(theta > 0)){
     0
   }else{
     -Inf
   }
 },
 logLikeApproxEval = function(y, thetaProposal, uProposal){
   if( (thetaProposal[1] > thetaProposal[2]) & (all(thetaProposal > 0))){
     xSamps <- uProposal*sqrt(thetaProposal[2])
     logCondLikes <- sapply(xSamps,
                           function(xsamp) {
                             sum(dnorm(y,
                                       xsamp,
                                       sqrt(thetaProposal[1] - thetaProposal[2]),
                                       log = TRUE)) })
     m <- max(logCondLikes)
     log(sum(exp(logCondLikes - m))) + m - log(length(y))
   }else{
     -Inf
   }
 },
 realY, numImportanceSamps, numMCMCIters, .99, recordEveryTh
)
res <- sampler(realParams)