Commit 87f02d7a authored by Chris Jewell's avatar Chris Jewell
Browse files

Deleted superfluous (and misleading) test files.

parent 84cb3605
#######################################
# Name: test.R #
# Created: 2016-06-14 #
# Author: Chris Jewell #
# Copyright: Chris Jewell 2016 #
# Purpose: Package testing functions #
# License: GPLv3 #
#######################################
riLogCP <- function(y,k,alpha,r,q) {
lambda <- q * r %*% (alpha * k)
sum(dpois(y,lambda, log = T))
}
testcampy <- function()
{
if (!require(dplyr, quietly = T))
stop("'test' requires package 'dplyr'")
data(campy)
campy <- filter(campy, Time == 1 & Location == 'A') %>%
filter(ChickenA + ChickenB + ChickenC + Bovine + Ovine + Environment > 0)
y <- campy$Human
y <-
array(
y, dim = c(length(y),1,1), dimnames = list(
type = campy$Type, time = '1', location = 'A'
)
)
X <-
campy[,c("ChickenA","ChickenB","ChickenC","Bovine","Ovine","Environment")]
for (i in 1:ncol(X))
X[,i] <- X[,i]# / sum(X[,i])
X <- array(
as.matrix(X), dim = c(nrow(X), ncol(X), 1),
dimnames = list(
type = campy$Type, source = c(
"ChickenA","ChickenB","ChickenC","Bovine","Ovine","Environment"
), time = '1'
)
)
k <-
array(
rep(1,dim(X)[2]), dim = c(dim(X)[2],1), dimnames = list(source = dimnames(X)$source, time = '1')
)
a_0 <- array(
1, dim = c(nrow(X), ncol(X), 1),
dimnames = list(
type = campy$Type, source = c(
"ChickenA","ChickenB","ChickenC","Bovine","Ovine","Environment"
), time = '1'
)
)
a_alpha <- rep(1, dim(X)[2])
list(y = y, R = X, time_names = "1", location_names = "A", source_names = colnames(X),
type_names = rownames(X),
k = k, a_q = 1, a_theta = 1, b_theta = 1, S_i = rep(1, dim(y)[1]), theta = 1, a_0 = a_0, a_alpha = a_alpha)
}
test_mcmc <- function(model, n.iter) {
nTypes <- length(model$y[[1]][[1]]$getData())
nSources <- length(model$alpha[[1]][[1]]$getData())
nTimes <- length(model$r)
nLocations <- length(model$alpha[[1]])
namesType <- names(model$y[[1]][[1]]$getData())
namesSource <- names(model$alpha[[1]][[1]]$getData())
namesTime <- names(model$y)
namesLocation <- names(model$y[[1]])
# Create posterior
posterior <- list()
posterior$alpha <- array(
dim = c(nSources,
nTimes,
nLocations,
n.iter),
dimnames = list(
source = namesSource, time = namesTime, location = namesLocation, iter = 1:n.iter
)
)
posterior$q <- array(
dim = c(
nTypes,
n.iter
),
dimnames = list(
type = namesType, iter = 1:n.iter
)
)
posterior$s <- array(
dim = c(
nTypes,
n.iter
),
dimnames = list(
type = namesType, iter = 1:n.iter
)
)
posterior$r <- array(
dim = c(nTypes,
nSources,
nTimes,
n.iter),
dimnames = list(
type = namesType, source = namesSource, time = namesTime, iter = 1:n.iter
)
)
# Create updaters for each time, location and source
r <- list()
alpha <- list()
for (time in 1:nTimes) {
alpha[[time]] <- list()
for (location in 1:nLocations) {
alpha[[time]][[location]] <- AdaptiveDirMRW$new(model$alpha[[time]][[location]], tune = c(0.4,0.4,0.4,0.4,0.4,0.1)) # User specificed starting variance?
}
r[[time]] <- list()
for (source in 1:nSources) { # TODO: change to an lapply?
r[[time]][[source]] <- AdaptiveLogDirMRW$new(model$r[[time]][[source]], toupdate = function() sample(nTypes, nTypes %/% 10), tune = rep(0.01, nTypes))
}
}
q <- PoisGammaDPUpdate$new(model$q)
# Save values to posterior
save_chain_state <- function(iter, posterior) {
for (time in 1:nTimes) {
for (location in 1:nLocations) {
posterior$alpha[, time, location, iter] <- model$alpha[[time]][[location]]$getData()
}
for (source in 1:nSources) { # TODO: change to an lapply?
posterior$r[ , source, time, iter] <- model$r[[time]][[source]]$getData()
}
}
posterior$q[, iter] <- model$q$getData()
posterior$s[, iter] <- model$q$s
posterior
}
posterior <- save_chain_state(1, posterior)
# Bit of user feedback and main MCMC loop
pb <- txtProgressBar(max = n.iter, style = 3)
for (i in 2:n.iter) {
# Do the updates
for (time in 1:nTimes) {
for (location in 1:nLocations) {
alpha[[time]][[location]]$update()
}
for (source in 1:nSources) { # TODO: change to an lapply?
r[[time]][[source]]$update()
}
}
q$update()
# Save chain state (possibly use a function?)
posterior <- save_chain_state(i, posterior)
setTxtProgressBar(pb, i)
}
#message("Alpha acceptance:")
#print(alpha$acceptance())
message("Summary R acceptance:")
for (time in 1:nTimes) {
acceptance_r_summary <- sapply(r[[time]], function(x) summary(x$acceptance()))
colnames(acceptance_r_summary) <- namesSource
print(acceptance_r_summary)
}
#print(sapply(r, function(x) x$acceptance()))
posterior
}
algtest <- function(n.iter = 1000)
{
data <- testcampy()
# y, R, Time, Location, Source, Type, Prev, a_q, a_theta, b_theta, a_0, a_alpha, S_i, theta
model <- DPModel(data$y, data$R,
data$time_names, data$location_names, data$source_names, data$type_names,
data$k, data$a_q, data$a_theta, data$b_theta, data$a_0,
data$a_alpha, data$S_i, data$theta)
rv <- NULL
t <- system.time(rv <- test_mcmc(model, n.iter))
message('Timing results')
print(t)
list(post = rv, model = model)
}
# res <- algtest(100)
# X11()
# par(mfrow = c(2, 4))
# nums <- cbind(sample(1 : dim(res$post$r)[1], 8), sample(1 : dim(res$post$r)[2], 8, replace = T))
# lapply(1 : 8, function(x) {
# plot(res$post$r[nums[x, 1], nums[x, 2], 1, ], type = "l")
# })
#
#
# X11()
# par(mfrow = c(2, 4))
# nums <- sample(1 : dim(res$post$q)[1], 8)
# lapply(nums, function(x) {
# plot(res$post$q[x, ], type = "l")
# })
#
# X11()
# par(mfrow = c(2, 3))
# lapply(1 : 6, function(x) {
# plot(res$post$alpha[x, 1, 1, ], type = "l", main = dimnames(res$post$alpha)[[1]][x])
# })
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