Title: | Constructs a Correlated Pseudo-Marginal Sampler |
---|---|
Description: | The primary function makeCPMSampler() generates a sampler function which performs the correlated pseudo-marginal method of Deligiannidis, Doucet and Pitt (2017) <arXiv:1511.04992>. If the 'rho=' argument of makeCPMSampler() is set to 0, then the generated sampler function performs the original pseudo-marginal method of Andrieu and Roberts (2009) <DOI:10.1214/07-AOS574>. The sampler function is constructed with the user's choice of prior, parameter proposal distribution, and the likelihood approximation scheme. Note that this algorithm is not automatically tuned--each one of these arguments must be carefully chosen. |
Authors: | Taylor Brown [aut, cre] |
Maintainer: | Taylor Brown <[email protected]> |
License: | MIT + file LICENSE |
Version: | 1.0.1 |
Built: | 2025-01-19 04:45:42 UTC |
Source: | https://github.com/cran/cPseudoMaRg |
checks if a log-density evaluation is not a valid number
isBadNum(num)
isBadNum(num)
num |
evaluation of a log-density |
TRUE or FALSE
isBadNum(NaN)
isBadNum(NaN)
correlated pseudo-marginal: generates functions that output a big vector
makeCPMSampler( paramKernSamp, logParamKernEval, logPriorEval, logLikeApproxEval, yData, numU, numIters, rho = 0.99, storeEvery = 1, nansInLLFatal = TRUE )
makeCPMSampler( paramKernSamp, logParamKernEval, logPriorEval, logLikeApproxEval, yData, numU, numIters, rho = 0.99, storeEvery = 1, nansInLLFatal = TRUE )
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 |
vector of theta samples
# 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)
# 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)
calculates the posterior mean point estimate
## S3 method for class 'cpmResults' mean(x, ...)
## S3 method for class 'cpmResults' mean(x, ...)
x |
a cpmResults object |
... |
arguments to be passed to or from methods. |
a vector of parameter estimates (posterior mean)
plots a cpmResults object
## S3 method for class 'cpmResults' plot(x, ...)
## S3 method for class 'cpmResults' plot(x, ...)
x |
a cpmResults object |
... |
arguments to be passed to or from methods. |
prints a cpmResults object
## S3 method for class 'cpmResults' print(x, ...)
## S3 method for class 'cpmResults' print(x, ...)
x |
a cpmResults object |
... |
arguments to be passed to or from methods. |
the same cpmResults object