inst/doc/introduction.R In mcunit: Unit Tests for MC Methods

```## -----------------------------------------------------------------------------
require(mcunit)
set.seed(10)

## -----------------------------------------------------------------------------
gibbsUpdate <- function(y, theta_j) {
# Samples theta_i given y and theta_j
mean <- 100 * (y - theta_j) / (100 + 0.1)
var <- 1. / (1. / 100 + 1. / 0.1)
rnorm(1, mean=mean, sd=sqrt(var))
}

## -----------------------------------------------------------------------------
randomScan <- function(theta, y, thinning) {
# Random Scan Gibbs
for(i in 1:thinning) {
# select index to update
i <- sample.int(2,1)
theta[i] <- gibbsUpdate(y, theta[i %% 2 + 1])
}
theta
}

## -----------------------------------------------------------------------------
obj <- list()
obj\$genprior <- function() rnorm(n=2, mean=0, sd=10)
obj\$gendata <- function(theta) sum(theta) + rnorm(n=1, mean=0, sd=sqrt(0.1))

## -----------------------------------------------------------------------------
priorDensity <- function(theta) prod(dnorm(theta, mean=0, sd=10))
likelihood <- function(theta, y) dnorm(y, mean=sum(theta), sd=sqrt(0.1))
testVec <- function(theta, y) c(theta[1], theta[2], priorDensity(theta), likelihood(theta, y))
obj\$test <- testVec

## -----------------------------------------------------------------------------
obj\$stepMCMC <- function(theta, dat, thinning) randomScan(theta, dat, thinning)
print(obj)

## -----------------------------------------------------------------------------
expect_mcmc(obj, thinning = 100)

## -----------------------------------------------------------------------------
expect_mcmc_reversible(obj, thinning = 10, nsteps = 10)

## -----------------------------------------------------------------------------
systematicScan <- function(theta, y, thinning) {
# Systematic Scan Gibbs
for(i in 1:thinning) {
theta[1] <- gibbsUpdate(y, theta[2])
theta[2] <- gibbsUpdate(y, theta[1])
}
theta
}

## -----------------------------------------------------------------------------
obj\$stepMCMC <- function(theta, dat, thinning) systematicScan(theta, dat, thinning)
expect_mcmc(obj, thinning = 100)

## -----------------------------------------------------------------------------
gibbsUpdate <- function(y, theta_j) {
# Samples theta_i given y and theta_j
mean <- 100 * (y - theta_j) / (100 + 0.1)
var <- 1. / (1. / 10 + 1. / sqrt(0.1))
rnorm(1, mean=mean, sd=sqrt(var))
}

## -----------------------------------------------------------------------------
obj\$stepMCMC <- function(theta, dat, thinning) randomScan(theta, dat, thinning)
```

Try the mcunit package in your browser

Any scripts or data that you put into this service are public.

mcunit documentation built on April 2, 2021, 5:06 p.m.