Simpler Nimble Dai Model

A simpler Nimble model

library("nimble")
library("regimeshifts")
library("tidyr")
library("ggplot2")
library("dplyr")

LSN version

A further simplified version:

• fix sigma_x to unity and theta to zero, reflecting the detrending and scaling.
• constrain sigma_y with a tight prior

code <- nimbleCode({
# sigma_x ~ dunif(1e-10, 1e2)
  
  ## highly constrained
  sigma_y ~ dgamma(10, 1000)

  ## Uninformative
       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] * x[i]
    x[i+1] ~ dnorm(mu_x[i], sd = 1) 
    y[i+1] ~ dnorm(y[i] + m * t[i] / t[N], sd = sigma_y) 
    
  }
})

Generate the test data:

set.seed(1000)
N <- 50
DF <- seq(0, 2000, length=N) # schedule for env degredation (increased dilution)
x <- numeric(N)   
x[1] <- 1.76e5 # initial density
   
for(day in 1:(N-1)){
 x[day+1] <- dai(x[day], DF = DF[day])
}

raw <- data.frame(t = 1:N, x = x)

Detrend:

raw$t <- 1:N
detrend <- loess(x ~ t, raw)
sigma <- sqrt(var(detrend$residuals))
data <- data.frame(x = detrend$residuals/sigma)
qplot(raw$t, data$x, geom='line')

constants <- list(N = N, t = raw$t)
inits <- list(m = 0, sigma_y = .01, y = rep(1,N))
thin <- 1e2
n_iter <- 1e6

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)

NULL

samples <- as.data.frame(as.matrix(Cmcmc$mvSamples))
samples <- samples[,1:(length(inits) - 1)]
df <- gather(samples)

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       m -0.02118436 0.02327253
2 sigma_y  0.01001409 0.00316794

Traces

ggplot(sample_n(df, 2e2)) + 
  geom_line(aes(seq_along(value), value)) + 
  facet_wrap(~key, scale='free')

Posteriors

ggplot(df) + 
  geom_density(aes(value)) + 
  facet_wrap(~key, scale='free')

Block sampler

constants <- list(N = N, t = raw$t)
inits <- list(m = 0, sigma_y = .01, y = rep(1,N))
Rmodel <- nimbleModel(code = code, constants = constants, data = data, inits = inits)
Cmodel <- compileNimble(Rmodel)
mcmcspec <- configureMCMC(Rmodel, print=FALSE,thin=thin)
mcmcspec$addSampler("RW_block", list(targetNodes=c('m','sigma_y'), adaptInterval=100))

[53] RW_block sampler;   targetNodes: m, sigma_y,  adaptive: TRUE,  adaptScaleOnly: FALSE,  adaptInterval: 100,  scale: 1,  propCov: identity

Rmcmc <- buildMCMC(mcmcspec)
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)
Cmcmc$run(n_iter)

NULL

samples <- as.data.frame(as.matrix(Cmcmc$mvSamples))
samples <- samples[,1:(length(inits)-1)]
df <- gather(samples)

df <- my_mcmc(code=lsn, constants, data, inits)

Error in eval(expr, envir, enclos): could not find function "my_mcmc"

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       m -0.025367634 0.027055911
2 sigma_y  0.009965489 0.003102902

Traces

ggplot(sample_n(df, 1e3)) + 
  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.3 RC (2015-03-06 r67947)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Debian GNU/Linux 8 (jessie)

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-1            yaml_2.1.13            
[7] knitr_1.9              

loaded via a namespace (and not attached):
 [1] assertthat_0.1   codetools_0.2-11 colorspace_1.2-6
 [4] DBI_0.3.1        digest_0.6.8     evaluate_0.5.5  
 [7] formatR_1.0      grid_3.1.3       gtable_0.1.2    
[10] igraph_0.7.1     labeling_0.3     lazyeval_0.1.10 
[13] magrittr_1.5     MASS_7.3-39      munsell_0.4.2   
[16] parallel_3.1.3   plyr_1.8.1       proto_0.3-10    
[19] Rcpp_0.11.5      reshape2_1.4.1   scales_0.2.4    
[22] stringr_0.6.2    tools_3.1.3