Nimble On Dai Model

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.01975407

Traces

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.0410488393 Traces 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.045674221

Traces

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
[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