In [15]:

    

library(nimble)


    

In [16]:

    

myBUGScode <- nimbleCode({mu ~ dnorm(0, sd = 100) 
                          sigma ~ dunif(0, 100) 
                          for(i in 1:10) y[i] ~ dnorm(mu, sd = sigma)})


    

In [17]:

    

myModel <- nimbleModel(myBUGScode)


    

    




defining model...
building model...
running calculate on model (any error reports that follow may simply reflect missing values in model variables) ... 
checking model sizes and dimensions... note that missing values (NAs) or non-finite values were found in model variables: mu, sigma, y. This is not an error, but some or all variables may need to be initialized for certain algorithms to operate properly.
model building finished.

In [18]:

    

library(igraph)
plot(myModel$getGraph())


    

In [19]:

    

myData <- rnorm(10, mean = 12, sd = 2)


    

In [20]:

    

myModel$setData(list(y = myData))


    

In [21]:

    

myModel$setInits(list(mu = 0, sigma = 1))


    

In [22]:

    

myMCMC <- buildMCMC(myModel)


    

In [23]:

    

compiled <- compileNimble(myModel, myMCMC)


    

    




compiling... this may take a minute. Use 'showCompilerOutput = TRUE' to see C++ compiler details.
compilation finished.

In [24]:

    

compiled$myMCMC$run(10000)


    

    




|-------------|-------------|-------------|-------------|
|-------------------------------------------------------|






    






NULL

In [10]:

    

samples <- as.matrix(compiled$myMCMC$mvSamples)
plot(density(samples[,'mu']))
plot(density(samples[,'sigma']))


    

In [40]:

    

plot(samples[ , 'mu'], type = 'l', xlab = 'iteration',  ylab = expression(mu))
plot(samples[ , 'sigma'], type = 'l', xlab = 'iteration', ylab = expression(sigma))


    

In [41]:

    

## Figure showing ACFs
acf1 <- with(acf(samples[, 'mu'], plot = FALSE), data.frame(lag, acf)) ## plot autocorrelation of mu sample
acf1$param <- 'mu'
acf2 <- with(acf(samples[, 'sigma'], plot = FALSE), data.frame(lag, acf))
acf2$param <- 'sigma'

acfall <- do.call('rbind', list(acf1, acf2))
acfall <- acfall[acfall$lag <= 20, ]

## line plot version
acfplot <- ggplot(data=acfall, mapping=aes(x=lag, y=acf)) +
    geom_line(stat = "identity", position = "identity") + facet_wrap(~param, ncol = 1)  + theme(legend.position="top", legend.title = element_blank(), text = element_text(size = 8)) + ylab('MCMC autocorrelation')               
print(acfplot)


    

In [36]:

    

pdf("/Users/blaine-mooers/Documents/Databases/LaTexWorkbench/BeamerSlideshows/2016Dec6orNov1_SCUG/mu.pdf",width = 4, height = 3.5)
plot(density(samples[,'mu']))
dev.off()


    

    





pdf: 2

In [37]:

    

pdf("/Users/blaine-mooers/Documents/Databases/LaTexWorkbench/BeamerSlideshows/2016Dec6orNov1_SCUG/sigma.pdf",width = 4, height = 3.5)
plot(density(samples[,'sigma']))
dev.off()


    

    





pdf: 2

In [26]:

    

library(ggplot2);library(ggExtra)


    

In [38]:

    

df1 <- data.frame(mu=samples[,'mu'], sigma=samples[,'sigma'])
(p1 <- ggplot(df1, aes(mu, sigma)) + geom_point() + theme_bw())
ggMarginal(p1)
ggMarginal(p1, type = "histogram", xparams = list(binwidth = 1, fill = "orange"))