Skip to content

Commit

Permalink
Merge branch 'master' of github.com:cboettig/labnotebook
Browse files Browse the repository at this point in the history
  • Loading branch information
cboettig committed Jan 1, 2015
2 parents d027f6e + 7fcee6e commit c5c20af
Show file tree
Hide file tree
Showing 33 changed files with 1,871 additions and 1,562 deletions.
1 change: 0 additions & 1 deletion .gitignore
Original file line number Diff line number Diff line change
Expand Up @@ -3,4 +3,3 @@ _gh-pages
google_analytics.json
.Rhistory
_cache
assets/figures/posts
3 changes: 3 additions & 0 deletions Makefile
Original file line number Diff line number Diff line change
@@ -1,6 +1,9 @@
## be sure to do `source ../.notebook-env.sh` fist
## Cannot run this in makefile since it doesn't understand 'source'
## Cannot run with bash -c "source ../.notebook-env.sh" since then it's not available to the parent shell
all:
make build
make deploy

build:
docker run --name labnotebook \
Expand Down
6 changes: 3 additions & 3 deletions _config.yml
Original file line number Diff line number Diff line change
Expand Up @@ -48,7 +48,7 @@ exclude:
- Dockerfile
- build.R
- _cache
- assets/figures/posts
knitr:
dev: svg
- Makefile
- deploy.sh
- fig.yml

2 changes: 0 additions & 2 deletions _posts/2014-12-03-nimble-explore.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,6 @@ library(methods)
library(knitr)
# set dev to png for graphs with > 1e3 points
opts_chunk$set(dev='png', fig.width=3, fig.height=4, cache=TRUE)
opts_knit$set(upload.fun = image_uri)
```

Expand Down Expand Up @@ -145,4 +144,3 @@ ylab = expression(beta))
Note the poor mixing (which is improved by the block sampler, as shown in the manual).


------
164 changes: 141 additions & 23 deletions _posts/2014-12-03-nimble-explore.md
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,10 @@ Working through quick-start example in [nimble manual](http://r-nimble.org/manua

The manual gives essentially no introduction to what appears to be a classic BUGS example model for stochastically failing pumps.

library(nimble)pumpCode <- modelCode({

```r
library(nimble)
pumpCode <- modelCode({
for (i in 1:N){
theta[i] ~ dgamma(alpha,beta)
lambda[i] <- theta[i]*t[i]
Expand All @@ -22,21 +25,37 @@ library(nimble)pumpCode <- modelCode({
alpha ~ dexp(1.0)
beta ~ dgamma(0.1,1.0)
})
```


```r
pumpConsts <- list(N = 10,
t = c(94.3, 15.7, 62.9, 126, 5.24,
31.4, 1.05, 1.05, 2.1, 10.5))pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))
31.4, 1.05, 1.05, 2.1, 10.5))
pumpData <- list(x = c(5, 1, 5, 14, 3, 19, 1, 1, 4, 22))
```


```r
pumpInits <- list(alpha = 1,
beta = 1,
theta = rep(0.1, pumpConsts$N))
```



```r
pump <- nimbleModel(code = pumpCode,
name = 'pump',
constants = pumpConsts,
data = pumpData,
inits = pumpInits)pump$getNodeNames() [1] "alpha" "beta" "lifted_d1_over_beta"
inits = pumpInits)

pump$getNodeNames()
```

```
[1] "alpha" "beta" "lifted_d1_over_beta"
[4] "theta[1]" "theta[2]" "theta[3]"
[7] "theta[4]" "theta[5]" "theta[6]"
[10] "theta[7]" "theta[8]" "theta[9]"
Expand All @@ -47,70 +66,169 @@ pump <- nimbleModel(code = pumpCode,
[25] "x[2]" "x[3]" "x[4]"
[28] "x[5]" "x[6]" "x[7]"
[31] "x[8]" "x[9]" "x[10]"

```

Note that we can see `theta` has our initial conditions, while `lambda` has not yet been initialized:

pump$theta [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
pump$lambda [1] NA NA NA NA NA NA NA NA NA NA

```r
pump$theta
```

```
[1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
```

```r
pump$lambda
```

```
[1] NA NA NA NA NA NA NA NA NA NA
```

Hmm, initially we cannot simulate `theta` values though (or rather, we just get NaNs and warnings if we do). At the moment I'm not clear on why, though seems to be due to the lifted node:

simulate(pump, 'theta')pump$theta [1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN

```r
simulate(pump, 'theta')
pump$theta
```

```
[1] NaN NaN NaN NaN NaN NaN NaN NaN NaN NaN
```

pump$lifted_d1_over_beta[1] NA

```r
pump$lifted_d1_over_beta
```

```
[1] NA
```



If we calculate the log probability density of the determinstic dependencies of alpha and beta nodes (i.e. the lifted node) then we're okay:

set.seed(0) ## This makes the simulations here reproduciblecalculate(pump, pump$getDependencies(c('alpha', 'beta'), determOnly = TRUE))[1] 0

```r
set.seed(0) ## This makes the simulations here reproducible
calculate(pump, pump$getDependencies(c('alpha', 'beta'), determOnly = TRUE))
```

```
[1] 0
```

```r
simulate(pump, 'theta')
```

pump$theta [1] 1.79180692 0.29592523 0.08369014 0.83617765 1.22254365 1.15835525
[7] 0.99001994 0.30737332 0.09461909 0.15720154

```r
pump$theta
```

```
[1] 1.79180692 0.29592523 0.08369014 0.83617765 1.22254365 1.15835525
[7] 0.99001994 0.30737332 0.09461909 0.15720154
```

We still need to initialize lambda, e.g. by calculating the probability density on those nodes:

calculate(pump, 'lambda')[1] 0
pump$lambda [1] 168.9673926 4.6460261 5.2641096 105.3583839 6.4061287
[6] 36.3723548 1.0395209 0.3227420 0.1987001 1.6506161

```r
calculate(pump, 'lambda')
```

```
[1] 0
```

```r
pump$lambda
```

```
[1] 168.9673926 4.6460261 5.2641096 105.3583839 6.4061287
[6] 36.3723548 1.0395209 0.3227420 0.1987001 1.6506161
```

though not entirely clear to me why the guide prefers to do this as the dependencies of theta (which clearly include lambda, but also other things). Also not clear if these `calculate` steps are necessary to proceed with the `MCMCspec` and `buildMCMC`, or compile steps. Let's reset the model[^1] and find out:

[^1]: Not completely certain that this destroys anything connected to the object as C pointers from before, but seems like it should.


```r
pump <- nimbleModel(code = pumpCode,
name = 'pump',
constants = pumpConsts,
data = pumpData,
inits = pumpInits)pump$theta [1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
pump$lambda [1] NA NA NA NA NA NA NA NA NA NA
inits = pumpInits)

pump$theta
```

```
[1] 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
```

```r
pump$lambda
```

```
[1] NA NA NA NA NA NA NA NA NA NA
```


Good, we're reset. Now we try:

Cpump <- compileNimble(pump)pumpSpec <- MCMCspec(pump)pumpSpec$addMonitors(c('alpha', 'beta', 'theta'))thin = 1: alpha, beta, theta

```r
Cpump <- compileNimble(pump)
pumpSpec <- MCMCspec(pump)
pumpSpec$addMonitors(c('alpha', 'beta', 'theta'))
```

```
thin = 1: alpha, beta, theta
```



```r
pumpMCMC <- buildMCMC(pumpSpec)
CpumpMCMC <- compileNimble(pumpMCMC, project = pump)
```

pumpMCMC <- buildMCMC(pumpSpec)CpumpMCMC <- compileNimble(pumpMCMC, project = pump)

CpumpMCMC(1000)NULL
```r
CpumpMCMC(1000)
```

```
NULL
```

```r
samples <- as.matrix(nfVar(CpumpMCMC, 'mvSamples'))
```


```r
plot(samples[ , 'alpha'], type = 'l', xlab = 'iteration',
ylab = expression(alpha))plot(samples[ , 'beta'], type = 'l', xlab = 'iteration',
ylab = expression(beta))plot(samples[ , 'alpha'], samples[ , 'beta'], xlab = expression(alpha),
ylab = expression(beta))assets/figures/posts/2014-12-03-nimble-explore/unnamed-chunk-15-1.pngassets/figures/posts/2014-12-03-nimble-explore/unnamed-chunk-15-2.pngassets/figures/posts/2014-12-03-nimble-explore/unnamed-chunk-15-3.png
ylab = expression(alpha))
plot(samples[ , 'beta'], type = 'l', xlab = 'iteration',
ylab = expression(beta))
plot(samples[ , 'alpha'], samples[ , 'beta'], xlab = expression(alpha),
ylab = expression(beta))
```

![plot of chunk unnamed-chunk-15](/assets/figures/posts/2014-12-03-nimble-explore/unnamed-chunk-15-1.png) ![plot of chunk unnamed-chunk-15](/assets/figures/posts/2014-12-03-nimble-explore/unnamed-chunk-15-2.png) ![plot of chunk unnamed-chunk-15](/assets/figures/posts/2014-12-03-nimble-explore/unnamed-chunk-15-3.png)

Note the poor mixing (which is improved by the block sampler, as shown in the manual).


------
3 changes: 2 additions & 1 deletion _posts/2014-12-03-ou-model-in-nimble.Rmd
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,8 @@ tags:
Sanity test with a simple model, Start with some sample data from an OU process:

```{r}
library(sde)
library("sde")
library("nimble")
set.seed(123)
d <- expression(0.5 * (10-x))
s <- expression(1)
Expand Down
Loading

0 comments on commit c5c20af

Please sign in to comment.