Testing nimble method on Dai model
(long data series)
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 <- 500
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.2075146 0.02739447
2 sigma 0.6156283 0.01975407Traces
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.003984249 0.0005859252
2 sigma_x 0.575993873 0.0185972953
3 theta -0.010341240 0.0410488393Traces
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.003001092 0.004531153
2 sigma_x 0.563929151 0.020224310
3 sigma_y 0.048817421 0.036018516
4 theta 0.010243795 0.045674221Traces
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.2.1 (2015-06-18)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux stretch/sid
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.2 ggplot2_1.0.1
[3] tidyr_0.2.0 regimeshifts_0.0.0.9000
[5] nimble_0.4 yaml_2.1.13
[7] knitr_1.10.5
loaded via a namespace (and not attached):
[1] igraph_1.0.1 Rcpp_0.12.0 magrittr_1.5
[4] MASS_7.3-43 munsell_0.4.2 colorspace_1.2-6
[7] R6_2.1.0 stringr_1.0.0 plyr_1.8.3
[10] tools_3.2.1 parallel_3.2.1 grid_3.2.1
[13] gtable_0.1.2 DBI_0.3.1 lazyeval_0.1.10
[16] digest_0.6.8 assertthat_0.1 reshape2_1.4.1
[19] formatR_1.2 codetools_0.2-14 evaluate_0.7
[22] labeling_0.3 stringi_0.5-5 scales_0.2.5
[25] proto_0.3-10