Commit 9b1cd80c authored by Chris Jewell's avatar Chris Jewell
Browse files

Added test for Dirichlet prior

parent 5af67f5d
^.*\.Rproj$
^\.Rproj\.user$
R/testDirichlet.R
R/test.R
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