MyNixOS website logo
Description

Constructs a Correlated Pseudo-Marginal Sampler.

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.

cPseudoMaRg

R-CMD-check

An implementation of the Correlated Pseudo-Marginal Sampler.

Install

Install from CRAN by typing

install.packages("cPseudoMaRg")

in an R console. Alternatively, install from Github by typing

devtools::install_github("tbrown122387/cpm")

Example

Another Random Effects Model that mimics the example in the above paper. They estimate a mean parameter, whereas the unknown parameters here are variance parameters. Also, this model's likelihood is nonidentifiable.

# y | x, theta ~ Normal(x, SSy)
# x | theta ~ Normal(0, SSx)
# theta = (SSy + SSx, SS_x)
# p(theta | y) propto p(y | theta)p(theta)
# approx p(y | theta) with mean( p(y | xi, theta)  ) where xi ~ p(xi | theta)

# real data
realxVar <- .2
realyVar <- .3
realTheta1 <- realxVar + realyVar
realTheta2 <- realxVar
realParams <- c(realTheta1, realTheta2)
numObs <- 10
realX <- rnorm(numObs, mean = 0, sd = sqrt(realxVar))
realY <- rnorm(numObs, mean = realX, sd = sqrt(realyVar))

# tuning params
numImportanceSamps <- 1000
numMCMCIters <- 10000
randomWalkScale <- 1.5
recordEveryTh <- 1
myLLApproxEval <- 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 = T)) })
    m <- max(logCondLikes)
    log(sum(exp(logCondLikes - m))) + m - log(length(y))
  }else{
    -Inf
  }
}
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 = myLLApproxEval,
  realY, numImportanceSamps, numMCMCIters, .99, recordEveryTh
)
res <- sampler(realParams)

# look at output
print(res)
plot(res)
Metadata

Version

1.0.1

License

Unknown

Platforms (75)

    Darwin
    FreeBSD
    Genode
    GHCJS
    Linux
    MMIXware
    NetBSD
    none
    OpenBSD
    Redox
    Solaris
    WASI
    Windows
Show all
  • aarch64-darwin
  • aarch64-genode
  • aarch64-linux
  • aarch64-netbsd
  • aarch64-none
  • aarch64_be-none
  • arm-none
  • armv5tel-linux
  • armv6l-linux
  • armv6l-netbsd
  • armv6l-none
  • armv7a-darwin
  • armv7a-linux
  • armv7a-netbsd
  • armv7l-linux
  • armv7l-netbsd
  • avr-none
  • i686-cygwin
  • i686-darwin
  • i686-freebsd
  • i686-genode
  • i686-linux
  • i686-netbsd
  • i686-none
  • i686-openbsd
  • i686-windows
  • javascript-ghcjs
  • loongarch64-linux
  • m68k-linux
  • m68k-netbsd
  • m68k-none
  • microblaze-linux
  • microblaze-none
  • microblazeel-linux
  • microblazeel-none
  • mips-linux
  • mips-none
  • mips64-linux
  • mips64-none
  • mips64el-linux
  • mipsel-linux
  • mipsel-netbsd
  • mmix-mmixware
  • msp430-none
  • or1k-none
  • powerpc-netbsd
  • powerpc-none
  • powerpc64-linux
  • powerpc64le-linux
  • powerpcle-none
  • riscv32-linux
  • riscv32-netbsd
  • riscv32-none
  • riscv64-linux
  • riscv64-netbsd
  • riscv64-none
  • rx-none
  • s390-linux
  • s390-none
  • s390x-linux
  • s390x-none
  • vc4-none
  • wasm32-wasi
  • wasm64-wasi
  • x86_64-cygwin
  • x86_64-darwin
  • x86_64-freebsd
  • x86_64-genode
  • x86_64-linux
  • x86_64-netbsd
  • x86_64-none
  • x86_64-openbsd
  • x86_64-redox
  • x86_64-solaris
  • x86_64-windows