Commit 614b5ad1 authored by Poppy Miller's avatar Poppy Miller
Browse files

Merging my updates to branch 'interface' of...

Merging my updates to branch 'interface' of fhm-chicas-code.lancs.ac.uk:jewell/sourceR-G into interface.
parents 61a628a2 5b143239
^.*\.Rproj$
^\.Rproj\.user$
R/testDirichlet.R
R/test.R
......@@ -134,7 +134,7 @@
#' data(campy)
#' zero_rows <- which(apply(campy[,c(2:7)], 1, sum) == 0)
#' campy <- campy[-zero_rows,]
#' prevs <- data.frame(Value = c(239/181, 196/113, 127/109,
#' prevs <- data.frame(Value = 1/c(239/181, 196/113, 127/109,
#' 595/97, 552/165, 524/86),
#' Source = colnames(campy[, 2:7]),
#' Time = rep(1, 6),
......@@ -489,7 +489,8 @@ HaldDP <- R6::R6Class(
} else {
private$r_updaters[[time]][[sources]] <- AdaptiveLogDirMRW$new(private$DPModel_impl$rNodes[[time]][[sources]],
toupdate = function() sample(private$nTypes, private$n_r),
tune = rep(0.01, private$nTypes))
tune = rep(1.0, private$nTypes),
batchsize=10)
}
}
}
......
......@@ -7,6 +7,65 @@
# class DAG model. #
##################################
AdaptiveSingleSiteMRW <- R6::R6Class(
"AdaptiveSingleSiteMRW",
public = list(
toupdate = NA,
naccept = NA,
ncalls = NA,
acceptbatch = NA,
batchsize = NA,
node = NA,
tune = NA,
initialize = function(node, tune = 0.1,
batchsize = 50) {
self$naccept <- 0
self$ncalls <- 0
self$acceptbatch <- 0
self$batchsize <- batchsize
self$tune <- tune
self$node <- node
},
update = function() {
self$ncalls <- self$ncalls + 1
old_data <- self$node$getData()
picur <- self$node$logPosterior()
# Propose using MHRW
self$node$data <- old_data * exp( rnorm(1, 0, self$tune) )
pican <- self$node$logPosterior()
alpha <- pican - picur + log(self$node$data/old_data)
if (is.finite(alpha) &
log(runif(1)) < alpha) {
self$acceptbatch <- self$acceptbatch + 1
}
else {
self$node$data <- old_data
}
private$adapt()
},
acceptance = function() {
self$naccept / self$ncalls
}
),
private = list(
adapt = function() {
if ((self$ncalls > 0) &
(self$ncalls %% self$batchsize == 0)) {
m <- ifelse(self$acceptbatch / self$batchsize > 0.44, 1,-1)
self$tune <-
exp(log(self$tune) + m * min(0.05, 1.0 / sqrt(self$ncalls)))
if(!is.finite(self$tune)) self$tune <- 1e-8
self$naccept <-
self$naccept + self$acceptbatch
self$acceptbatch <- 0
}
}
)
)
#' AdaptiveMultiMRW
#'
#' This class implements an Adaptive Multi-site Metropolis random walk
......@@ -258,7 +317,7 @@ AdaptiveLogDirMRW <- R6::R6Class(
batchsize = 50) {
self$toupdate = toupdate
self$naccept <- rep(0, length(node$getData()))
self$ncalls <- 0
self$ncalls <- rep(0, length(node$getData()))
self$acceptbatch <-
rep(0, length(node$getData()))
self$batchsize <- batchsize
......@@ -266,8 +325,9 @@ AdaptiveLogDirMRW <- R6::R6Class(
self$node <- node
},
update = function() {
self$ncalls <- self$ncalls + 1
for (j in self$toupdate()) {
updIdx <- self$toupdate()
for (j in updIdx) {
self$ncalls[j] <- self$ncalls[j] + 1
old_data <- self$node$getData()
picur <- self$node$logPosterior()
......@@ -286,23 +346,94 @@ AdaptiveLogDirMRW <- R6::R6Class(
else {
self$node$data <- old_data
}
private$adapt(j)
}
if(self$ncalls %% 50 == 0) private$adapt()
},
acceptance = function() {
self$naccept <- self$naccept + self$acceptbatch
self$naccept / self$ncalls
}
),
private = list(
adapt = function() {
if ((self$ncalls > 0) &
(self$ncalls %% self$batchsize == 0)) {
m <- ifelse(self$acceptbatch / self$batchsize > 0.44, 1,-1)
self$tune <-
exp(log(self$tune) + m * min(0.05, 1.0 / sqrt(self$ncalls)))
self$naccept <-
self$naccept + self$acceptbatch
self$acceptbatch <- self$acceptbatch * 0
adapt = function(j) {
if ((self$ncalls[j] > 0) &
(self$ncalls[j] %% self$batchsize == 0)) {
m <- ifelse(self$acceptbatch[j] / self$batchsize > 0.44, 1,-1)
self$tune[j] <-
exp(log(self$tune[j]) + m * min(0.05, 1.0 / sqrt(self$ncalls[j])))
if(!is.finite(self$tune[j])) {
warning('AdaptiveLogDirMRW tuning variance non-finite')
self$tune[j] <- 1e-9
}
self$naccept[j] <-
self$naccept[j] + self$acceptbatch[j]
self$acceptbatch[j] <- 0
}
}
)
)
AdaptiveLogDirMRW2 <- R6::R6Class(
"AdaptiveDirMRW2",
public = list(
toupdate = NA,
naccept = NA,
ncalls = NA,
acceptbatch = NA,
batchsize = NA,
node = NA,
tune = NA,
initialize = function(node, toupdate = function() 1:length(node$data),
tune = rep(0.1, length(node$data)),
batchsize = 50) {
self$toupdate = toupdate
self$naccept <- rep(0, length(node$data))
self$ncalls <- rep(0, length(node$data))
self$acceptbatch <-
rep(0, length(node$data))
self$batchsize <- batchsize
self$tune <- tune
self$node <- node
},
update = function() {
updIdx <- self$toupdate()
for (j in updIdx) {
self$ncalls[j] <- self$ncalls[j] + 1
old_data <- self$node$data
picur <- self$node$logPosterior()
# Propose using MHRW
self$node$data[j] <- old_data[j] * exp( rnorm(1, 0, self$tune[j]) )
pican <- self$node$logPosterior()
alpha <- pican - picur + log(self$node$data[j]/old_data[j])
if (is.finite(alpha) &
log(runif(1)) < alpha) {
self$acceptbatch[j] <- self$acceptbatch[j] + 1
}
else {
self$node$data <- old_data
}
private$adapt(j)
}
},
acceptance = function() {
self$naccept <- self$naccept + self$acceptbatch
self$naccept / self$ncalls
}
),
private = list(
adapt = function(j) {
if ((self$ncalls[j] > 0) &
(self$ncalls[j] %% self$batchsize == 0)) {
m <- ifelse(self$acceptbatch[j] / self$batchsize > 0.44, 1,-1)
self$tune[j] <-
exp(log(self$tune[j]) + m * min(0.05, 1.0 / sqrt(self$ncalls[j])))
self$naccept[j] <-
self$naccept[j] + self$acceptbatch[j]
self$acceptbatch[j] <- 0
}
}
)
......
......@@ -71,7 +71,6 @@ DPModel_impl <- R6::R6Class(
# Construct y, lambda, a, R, and k for location/time pairs
for (time in unique(Time)) {
# Prevalences
k <- DataNode$new(data = setNames(prev[, time], Sources), name = paste('k',time,sep = '_'))
......
......@@ -208,7 +208,8 @@ PoissonNode <- R6::R6Class(
self$addParent(offset, 'offset')
},
logDensity = function() {
lambda <- sapply(self$parents, function(parent) parent$getData())
lambda <- as.data.frame(sapply(self$parents,
function(parent) parent$getData()))
lambda <- apply(lambda, 1, prod)
sum(dpois(self$data, lambda, log = T))
}
......@@ -290,6 +291,25 @@ DirichletNode <- R6::R6Class(
)
)
DirichletNode2 <- R6::R6Class(
"DirichletNode2",
inherit = StochasticNode,
public = list(
initialize = function(data, alpha, name) {
super$initialize(name = name)
self$data <- data
self$addParent(alpha, name='alpha')
},
logDensity = function() {
sum(dgamma(self$data, shape=self$parents$alpha$getData(), rate=1, log=T))
},
getData = function() {
self$data / sum(self$data)
}
)
)
#' Transformed Dirichlet node
#'
#' Uses transformation due to Betancourt 2013
......
testDirichlet <- function() {
data(campy)
X <- c(40, 10, 24, 1, 1, 5)
z <- c(2, 2, 4, 0, 0, 0)
prior <- DataNode$new(data = 1 + X,
'prior')
r <- X/sum(X)
rNode <- DirichletNode$new(data = r,
alpha = prior,
name='r')
a <- DataNode$new(data=0.1, name='a')
b <- DataNode$new(data=0.1, name='b')
beta <- GammaNode$new(data = 425,
shape = a,
rate = b,
name = 'beta')
y <- PoissonNode$new(data = z,
lambda = beta,
offset = rNode,
name = 'y')
# MCMC
posterior <- matrix(ncol = length(r) + 1, nrow = 1000)
posterior[1,1] <- beta$getData()
posterior[1,2:ncol(posterior)] <- rNode$getData()
rUpd <- AdaptiveLogDirMRW$new(node = rNode,
toupdate = function() sample(length(r), 3))
betaUpd <- AdaptiveSingleSiteMRW$new(node = beta)
for(i in 2:1000) {
betaUpd$update()
rUpd$update()
posterior[i,1] <- beta$getData()
posterior[i,2:ncol(posterior)] <- rNode$getData()
}
message("Beta accept:")
print(betaUpd$acceptance())
message("R accept: ")
print(rUpd$acceptance())
posterior
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment