test_bcplm.R
if (FALSE){
# load package
library(tweedie)
library(rbenchmark)
library(ggplot2)
library(coda)
library(lme4)
library(amer)
library(cplm)
options(error = recover)
#setwd("~/2011/cplm")
setwd("C:\\Documents and Settings\\CAB2007\\My Documents\\2012\\cplm")
load("./data/fineroot.RData")
source("./R/classMethods.R")
source("./R/cpglm.R")
source("./R/utilities.R")
#dyn.load("./src/cplm.so")
link="log"
control=list()
trace=T
n.chains=3
n.iter = 200
n.burnin = 0
n.sims=600
n.thin=1
n.report=2
inits=NULL
bound.p=c(1.2,1.9)
phi.shape=0.001
phi.scale=0.001
tune.iter=3000
n.tune= 15
tune.weight= 0.25
bound.phi =100
prior.beta.var=NULL
prior.beta.mean=NULL
contrasts = NULL
block.update.beta = FALSE
dyn.load("src/bcplm.dll")
.Call("init")
bcplm <- function(formula, link = "log", data, inits = NULL,
weights, offset, subset, na.action, contrasts = NULL,
n.chains = 3, n.iter = 2000, n.burnin = floor(n.iter / 2),
n.thin = max(1, floor(n.chains * (n.iter - n.burnin) / n.sims)),
n.sims = 1000, n.report = 1000, prior.beta.mean = NULL,
prior.beta.var = NULL, bound.phi = 100, bound.p = c(1.01, 1.99),
tune.iter = 4000, n.tune = 10, tune.weight = 0.25,
basisGenerators = c("tp","tpU","bsp","sp2d"),
block.update.beta = FALSE, ...) {}
set.seed(10)
# test glm
call <- amer:::expand.call(bcplm,call("bcplm", RLD~ Stock * Zone ,
data= FineRoot))
formula <- RLD ~ factor(Zone) * factor(Stock)
# test glmm
call <- amer:::expand.call(bcplm,call("bcplm",RLD ~ Stock * Zone + (1|Plant) ,
data = FineRoot))
formula <- RLD ~ Stock * Zone + (1|Plant)
fineroot$a <- rnorm(nrow(fineroot))
formula <- RLD~ Zone+(1+a|Plant)+(1|Stock)
contrasts = NULL
call <- amer:::expand.call(bcplm,call("bcplm", RLD~ Zone+(1+a|Plant)+(1|Stock),data=fineroot))
dyn.unload("src/bcplm.dll")
# fit the fineroot data with Bayesian models
# first use the latent variable approach
set.seed(10)
fit1 <- bcplm(RLD~ factor(Zone)*factor(Stock), data=fineroot,
n.iter=20000, n.burnin=10000,
n.thin=20, n.report=5000, method = "gibbs")
gelman.diag(fit1$sims.list)
# diagnostic plots
acfplot(fit1$sims.list,lag.max=20)
xyplot(fit1$sims.list)
densityplot(fit1$sims.list)
# summary results
summary(fit1$sims.list)
# now try direct density evaluation (This is much slower)
set.seed(10)
fit2 <- bcpglm(RLD~ factor(Zone)*factor(Stock), data=fineroot,
n.iter=11000, n.burnin=1000,
n.thin=10, n.report=5000)
gelman.diag(fit2$sims.list)
summary(fit2)
# now fit the Bayesian model to an insurance loss triangle
# (see Peters et al. 2009)
fit3 <- bcplm(increLoss~ factor(year)+factor(lag), data=InsReserve,
tune.iter=5000, n.tune = 10,
n.iter=11000, n.burnin=1000,
n.thin=10, n.report=5000, bound.p=c(1.1,1.95))
gelman.diag(fit3$sims.list)
summary(fit3$sims.list)
fit0<- cpglmm(RLD ~ Stock + (1|Plant), data=fineroot)
set.seed(10)
fit1 <- bcplm(RLD ~ Stock + (1|Plant), data=fineroot,
n.iter=60000, n.burnin=10000,
n.thin=50, n.report=5000, method = "gibbs")
gelman.diag(fit1$sims.list)
# diagnostic plots
acfplot(fit1$sims.list,lag.max=20)
xyplot(fit1$sims.list)
densityplot(fit1$sims.list)
# summary results
summary(fit1$sims.list)
set.seed(12)
fit2<- bcplm(RLD ~ Stock + (1|Plant), data=fineroot,
n.iter=20000, n.burnin=10000,
n.thin=20, n.report=5000, method = "gibbs")
summary(fit2$sims.list)
xyplot(fit2$sims.list)
mu1 <- as.numeric(exp(input@X %*% input@fixef + t(input@Zt) %*% input@u))
.Call("bcplm_update_mu", input)
sum(mu1 - input@mu)
#####################
# check the computation of full conditionals
# 1. post of beta_k
post_betak <- function(x, da){
beta <- da@fixef
p2 <- 2 - da@p
p1 <- da@p - 1
k <- da@k + 1
beta[k] <- x
mu <- as.numeric(exp(da@X %*% beta + t(da@Zt) %*% da@u))
-(sum(mu^p2)/p2 + sum(da@y * mu^(-p1))/p1)/da@phi -
0.5 * (x - da@pbeta.mean[k])^2/da@pbeta.var[k]
}
input@k <- as.integer(1)
post_betak(1.0, input)
.Call("bcplm_update_mu", input)
.Call("bcplm_post_betak", 1.0, input)
# 2. post of phi:
post_phi <- function(x, da){
sum(log(dtweedie(da@y, mu = da@mu, power = da@p, phi = x)))
}
post_phi(1.3, input)
.Call("bcplm_post_phi", 1.3, input)
# 2. post of p:
post_p <- function(x, da){
sum(log(dtweedie(da@y, mu = da@mu, power = x, phi = da@phi)))
}
post_p(1.3, input)
.Call("bcplm_post_p", 1.3, input)
# 4. post of u_k: assume "bcplm_update_mu" is called
post_uk <- function(x, da){
u <- da@u
p2 <- 2 - da@p
p1 <- da@p - 1
k <- da@k + 1
u[k] <- x
mu <- as.numeric(exp(da@X %*% da@fixef + t(da@Zt) %*% u))
-(sum(mu^p2)/p2 + sum(da@y * mu^(-p1))/p1)/da@phi -
0.5 * x^2/as.numeric(da@Sigma[[1]])
}
# uniform prior
post_sigma <- function(x, da){
dm <- da@dims
as.numeric(-dm["n.u"]/2 * log(x) - sum(da@u^2)/(2 * x^2))
}
post_sigma(0.08, input)
xx = seq(0.01, 0.2, by = 0.01)
yy <- sapply(xx, function(t) post_sigma(t, input))
yy2 <- yy - max(yy)
plot(xx, exp(yy2), type = "l")
input@k <- as.integer(1)
post_uk(1.0, input)
.Call("bcplm_update_mu", input)
.Call("bcplm_post_uk", 1.0, input)
post_phi(1.0, input)
.Call("bcplm_post_phi", 1.0, input)
pv <- c(0.01802352, 0.05255602, 0.04535603, 0.05694167, 0.11490249, 0.12897299,
0.002066605, 0.0006284337,
0.03479402, 0.0292, 0.0310, 0.0880, 0.02506620, 0.0257, 0.0347, 0.0387,
0.5)
input@mh.var <- pv
fit <- cpglmm(RLD ~ Stock * Zone + (1|Plant), data = FineRoot)
input@u <- fit@ranef
input@Sigma <- list(matrix(as.vector(VarCorr(fit)[[1]])))
do_mcmc <- function(input){
psd <- sqrt(input@mh.var)
ig.shape <- 0.0001
ig.scale <- 0.0001
dm <- input@dims
acc <- rep(0, dm["n.beta"] + dm["n.u"] + 2 + 1)
# update xb, zu, eta and mu slots etc
.Call("bcplm_update_mu", input)
sim <- matrix(0, dm["n.keep"], dm["n.beta"] + dm["n.u"] + 3)
ns <- 0
for (i in 1:dm["n.iter"]){
# sample beta
for (k in 1:dm["n.beta"]){
input@k <- as.integer(k - 1)
bk <- metrop_rw(1, input@fixef[k], psd[k], post_betak, input)
input@fixef[k] <- as.vector(bk)
acc[k] <- acc[k] + attr(bk, "accept")
}
# sample phi and p
input@mu <- as.numeric(exp(input@X %*% input@fixef + t(input@Zt) %*% input@u))
pos <- dm["n.beta"] + 1
phi <- metrop_rw(1, input@phi, psd[pos], post_phi, input, lower = 0, upper = 100)
input@phi <- phi
acc[pos] <- acc[pos] + attr(phi, "accept")
p <- metrop_rw(1, input@p, psd[pos + 1], post_p, input, lower = 1.01, upper = 1.99)
input@p <- p
acc[pos + 1] <- acc[pos + 1] + attr(p, "accept")
# sample u
for (k in 1:dm["n.u"]){
input@k <- as.integer(k - 1)
pos <- dm["n.beta"] + 2 + k
uk <- metrop_rw(1, input@u[k], psd[pos], post_uk, input)
input@u[k] <- as.vector(uk)
acc[pos] <- acc[pos] + attr(uk, "accept")
}
# simulate Sigma
#input@Sigma[[1]] <- matrix(1/rgamma(1, dm["n.u"]/2 + ig.shape, sum(input@u^2)/2 + ig.scale))
pos <- dm["n.beta"] + 2 + dm["n.u"] + 1
sigma <- metrop_rw(1, as.numeric(sqrt(input@Sigma[[1]])),
psd[pos], post_sigma, input, lower = 0, upper = 100)
input@Sigma <- list(matrix(sigma^2))
acc[pos] <- acc[pos] + attr(sigma, "accept")
# set parameters
if (i > dm["n.burnin"] && ((i - dm["n.burnin"])%%dm["n.thin"] == 0)){
ns <- ns + 1
sim[ns, ] <- c(input@fixef, input@phi, input@p, input@u, as.numeric(input@Sigma[[1]]))
}
}
return(list(accept = acc, sim = sim))
}
dm <- input@dims
dm["n.iter"] <- as.integer(500)
dm["n.burnin"] <- as.integer(0)
dm["n.thin"] <- as.integer(1)
dm["n.keep"] <- as.integer(500)
input@dims <- dm
dm <- input@dims
dm["n.iter"] <- as.integer(6000)
dm["n.burnin"] <- as.integer(1000)
dm["n.thin"] <- as.integer(5)
dm["n.keep"] <- as.integer(3000)
input@dims <- dm
a <- do_mcmc(input)
a[[1]]/500
out <- vector("list", dm["n.chains"])
for (i in 1:dm["n.chains"]){
init <- input@inits[[i]]
input@fixef <- init[1:dm["n.beta"]]
input@phi <- init[dm["n.beta"] + 1]
input@p <- init[dm["n.beta"] + 2]
input@u <- init[(dm["n.beta"] + 3):(dm["n.beta"] + dm["n.u"] + 2)]
input@Sigma[[1]] <- matrix(init[length(init)])
.Call("bcplm_update_mu", input)
cat("chain:")
cat(i)
out[[i]] <- do_mcmc(input)
}
a = lapply(out, "[[", 2)
b= lapply(a, as.mcmc)
b <- as.mcmc.list(b)
summary(b)
set.seed(10)
ss1 <- mh(500, postLog, xinit=5,lb=0,sigma=0.1)
ss2 <- mh(500, postLog, xinit=5,lb=0,sigma=50)
fun <- function(x, sd2 = 5){
dnorm(x, 0, sd = sd2, log = TRUE)
}
n <- 100000
set.seed(10)
s1 <- metrop_rw(n, par = 0, sd = 7, fun, sd2 = 5)
set.seed(10)
s2 <- metrop_rw1(n, par = 0, v = 50, fun, sd = 5)
c(mean(s1), sd(s1))
c(mean(s2), sd(s2))
attr(s1, "accept")/n
for (i in 1:20){
m1 <- MCMCmetrop1R(fun, 0, burnin = 0, seed =i,mcmc = n, V = matrix(50), sd = 5)
print(c(mean(m1), sd(m1)))
}
library(rbenchmark)
f1 <- function() metrop_rw1(n, par = 0, v = 50, fun, sd = 5)
f2 <- function() MCMCmetrop1R(fun, 0, burnin = 0, mcmc = n,
V = matrix(50), sd = 5, verbose = FALSE)
benchmark(f1(), f2(), replications = 50)
metrop_rw1(10, 0, 0.1, post_uk, input)
dyn.load("src/bcplm.dll")
.Call("init")
dyn.unload("src/bcplm.dll")
}
