Lsn Model Comparisons

library(knitr)
library(nimble)
library(sde)
library(ggplot2)
library(tidyr)
opts_chunk$set(dev = 'png', 
               fig.width = 5, 
               fig.height = 5, 
               results = 'hide')

Stick with SDE example as a reference point:

set.seed(123)
d <- expression(0.5 * (10-x))
s <- expression(1) 
data <- as.data.frame(sde.sim(X0=6,drift=d, sigma=s, T=20, N=100))

sigma.x not provided, attempting symbolic derivation.

ggplot(data) + geom_line(aes(seq_along(x),x))
N <- length(data$x)

LSN version

Modify the LSN model to explicitly model the changing parameter as a hidden, stochastic variable

lsn <- nimbleCode({
        theta ~ dunif(1e-10, 100.0)
        sigma_x ~ dunif(1e-10, 100.0)
        sigma_y ~ dunif(1e-10, 100.0)
                    m ~ dunif(-1e2, 1e2)
             x[1] ~ dunif(0, 100)
             y[1] ~ dunif(0, 100) 

  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]
    y[i+1] ~ dnorm(mu_y[i], sd = sigma_y) 
  }
})

Constants in the model definition are the length of the dataset, \(N\) and the time points of the sample. Note we’ve made time explicit, we’ll assume uniform spacing here.

constants <- list(N = N, t = 1:N)

Initial values for the parameters

inits <- list(theta = 6, m = 0, sigma_x = 1, sigma_y = 1, y = rep(1,N))

and here we go:

Rmodel <- nimbleModel(code = lsn, 
                      constants = constants, 
                      data = data, 
                      inits = inits)
Cmodel <- compileNimble(Rmodel)

mcmcspec <- configureMCMC(Rmodel, print=TRUE,thin=2e2)
Rmcmc <- buildMCMC(mcmcspec)
Cmcmc <- compileNimble(Rmcmc, project = Cmodel)

Cmcmc$run(1e6)

NULL

and examine results

samples <- as.data.frame(as.matrix(Cmcmc$mvSamples))
dim(samples)

[1] 5000  206

samples <- samples[,1:4]
long <- gather(samples)

apply(samples, 2, mean)

            m       sigma_x       sigma_y         theta 
-1.675333e-05  3.875174e-01  4.066872e-02  1.022867e+01 

ggplot(long) + 
    geom_line(aes(seq_along(value), value)) + 
    facet_wrap(~key, scale='free')

ggplot(long) + 
    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] stats4    splines   methods   stats     graphics  grDevices
[7] utils     datasets  base     

other attached packages:
 [1] tidyr_0.2.0   ggplot2_1.0.0 sde_2.0.13    zoo_1.7-11   
 [5] fda_2.4.4     Matrix_1.1-5  MASS_7.3-39   nimble_0.3-1 
 [9] yaml_2.1.13   knitr_1.9    

loaded via a namespace (and not attached):
 [1] codetools_0.2-11 colorspace_1.2-6 digest_0.6.8    
 [4] evaluate_0.5.5   formatR_1.0      grid_3.1.3      
 [7] gtable_0.1.2     igraph_0.7.1     labeling_0.3    
[10] lattice_0.20-30  munsell_0.4.2    plyr_1.8.1      
[13] proto_0.3-10     Rcpp_0.11.5      reshape2_1.4.1  
[16] scales_0.2.4     stringr_0.6.2    tools_3.1.3