Nothing
## -----------------------------------------------------------------------------
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)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.