Testing nimble method on Dai model
(shorter data series)
devtools::install_github("cboettig/regimeshifts")Define our mcmc procedure in Nimble
my_mcmc <- function(code, constants, data, inits, n_iter=1e5, thin = 1e2){
Rmodel <- nimbleModel(code = code, constants = constants, data = data, inits = inits)
Cmodel <- compileNimble(Rmodel)
mcmcspec <- configureMCMC(Rmodel, print=FALSE,thin=thin)
Rmcmc <- buildMCMC(mcmcspec)
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)
Cmcmc$run(n_iter)
samples <- as.data.frame(as.matrix(Cmcmc$mvSamples))
samples <- samples[,1:(length(inits)-1)]
gather(samples)
}Generate the test data:
set.seed(1000)
max_days <- 50
DF <- seq(0, 2000, length=max_days) # schedule for env degredation (increased dilution)
x <- numeric(max_days)
x[1] <- 1.76e5 # initial density
for(day in 1:(max_days-1)){
x[day+1] <- dai(x[day], DF = DF[day])
}
raw <- data.frame(t = 1:max_days, x = x)Detrend:
N <- length(raw$x)
raw$t <- 1:N
detrend <- loess(x ~ t, raw)
data <- data.frame(x = detrend$residuals/sqrt(var(detrend$residuals)))
qplot(raw$t, data$x, geom='line')
OU Model
ou <- nimbleCode({
theta ~ dunif(1e-10, 100.0)
r ~ dunif(1e-10, 20.0)
sigma ~ dunif(1e-10, 100)
x[1] ~ dunif(0, 100)
for(t in 1:(N-1)){
mu[t] <- x[t] + r * (theta - x[t])
x[t+1] ~ dnorm(mu[t], sd = sigma)
}
})
ou_constants <- list(N = N)
ou_inits <- list(theta = 0, r = 1e-3, sigma = 1)Run the mcmc
df <- my_mcmc(code=ou, ou_constants, data, ou_inits)Summary statistics
summarise(group_by(df, key), mean=mean(value), std=sqrt(var(value)))Source: local data frame [2 x 3]
key mean std
1 r 0.7866844 0.1600235
2 sigma 1.0375633 0.1113334Traces
ggplot(df) +
geom_line(aes(seq_along(value), value)) +
facet_wrap(~key, scale='free')
Posteriors
ggplot(df) +
geom_density(aes(value)) +
facet_wrap(~key, scale='free')
LSN version
A modified version of the LSN model to explicitly model the changing parameter as a hidden variable changing at constant rate
lsn <- nimbleCode({
theta ~ dunif(-1e2, 1e2)
sigma_x ~ dunif(1e-10, 1e2)
m ~ dunif(-1e2, 1e2)
x[1] ~ dunif(-1e3, 1e3)
y[1] ~ dunif(-1e3, 1e3)
for(i in 1:(N-1)){
mu_x[i] <- x[i] + y[i] * (theta - x[i])
x[i+1] ~ dnorm(mu_x[i], sd = sigma_x)
y[i+1] <- y[i] + m * t[i] / t[N]
}
})
constants <- list(N = N, t = raw$t)
inits <- list(theta = 0, m = 0, sigma_x = 1, y = rep(1,N))df <- my_mcmc(code=lsn, constants, data, inits)Summary statistics
summarise(group_by(df, key), mean=mean(value), std=sqrt(var(value)))Source: local data frame [3 x 3]
key mean std
1 m -0.02060616 0.02779506
2 sigma_x 1.03148569 0.11158787
3 theta -0.03136262 0.18565435Traces
ggplot(df) +
geom_line(aes(seq_along(value), value)) +
facet_wrap(~key, scale='free')
Posteriors
ggplot(df) +
geom_density(aes(value)) +
facet_wrap(~key, scale='free')
LSN, stochastic hidden variable
Define and model and run MCMC
lsn <- nimbleCode({
theta ~ dunif(-1e2, 1e2)
sigma_x ~ dunif(1e-10, 1e2)
sigma_y ~ dunif(1e-10, 1e2)
m ~ dunif(-1e2, 1e2)
x[1] ~ dunif(-1e3, 1e3)
y[1] ~ dunif(-1e3, 1e3)
for(i in 1:(N-1)){
mu_x[i] <- x[i] + y[i] * (theta - x[i])
x[i+1] ~ dnorm(mu_x[i], sd = sigma_x)
mu_y[i] <- y[i] + m * t[i] / t[N]
y[i+1] ~ dnorm(mu_y[i], sd = sigma_y)
}
})
constants <- list(N = N, t = raw$t)
inits <- list(theta = 0, m = 0, sigma_x = 1, sigma_y = 1, y = rep(1,N))df <- my_mcmc(code=lsn, constants, data, inits)Summary statistics
summarise(group_by(df, key), mean=mean(value), std=sqrt(var(value)))Source: local data frame [4 x 3]
key mean std
1 m -0.008170135 0.06970863
2 sigma_x 1.023650010 0.11068862
3 sigma_y 0.171747380 0.14767809
4 theta -0.025881301 0.18058435Traces
ggplot(df) +
geom_line(aes(seq_along(value), value)) +
facet_wrap(~key, scale='free')
Posteriors
ggplot(df) +
geom_density(aes(value)) +
facet_wrap(~key, scale='free')
sessionInfo()R version 3.1.2 (2014-10-31)
Platform: x86_64-pc-linux-gnu (64-bit)
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] methods stats graphics grDevices utils datasets
[7] base
other attached packages:
[1] dplyr_0.4.1 ggplot2_1.0.0
[3] tidyr_0.2.0 regimeshifts_0.0.0.9000
[5] nimble_0.3 yaml_2.1.13
[7] knitr_1.9
loaded via a namespace (and not attached):
[1] assertthat_0.1 bitops_1.0-6 codetools_0.2-10
[4] colorspace_1.2-4 DBI_0.3.1 devtools_1.7.0
[7] digest_0.6.8 evaluate_0.5.5 formatR_1.0
[10] grid_3.1.2 gtable_0.1.2 httr_0.6.1
[13] igraph_0.7.1 labeling_0.3 lazyeval_0.1.10
[16] magrittr_1.5 MASS_7.3-39 munsell_0.4.2
[19] parallel_3.1.2 plyr_1.8.1 proto_0.3-10
[22] Rcpp_0.11.4 RCurl_1.95-4.5 reshape2_1.4.1
[25] scales_0.2.4 stringr_0.6.2 tools_3.1.2