makeCPMSampler: correlated pseudo-marginal: generates functions that output a...

In tbrown122387/cpm: Constructs a Correlated Pseudo-Marginal Sampler

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)

tbrown122387/cpm documentation built on Nov. 22, 2022, 5:32 p.m.