Created
December 3, 2020 04:41
-
-
Save oliviergimenez/e41e9cb99174f2124f948308e19ca7ec to your computer and use it in GitHub Desktop.
Use Nimble to generate inits for latent states
We can make this file beautiful and searchable if this error is corrected: No commas found in this CSV file in line 0.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
h1;h2;h3;h4;h5;h6;h7;males;females | |
1;1;1;1;1;1;0;1;0 | |
1;1;1;1;1;0;0;0;1 | |
1;1;1;1;0;0;0;1;0 | |
1;1;1;1;0;0;0;0;1 | |
1;1;0;1;1;1;0;0;1 | |
1;1;0;0;0;0;0;1;0 | |
1;1;0;0;0;0;0;1;0 | |
1;1;0;0;0;0;0;1;0 | |
1;1;0;0;0;0;0;1;0 | |
1;1;0;0;0;0;0;0;1 | |
1;1;0;0;0;0;0;0;1 | |
1;0;1;0;0;0;0;1;0 | |
1;0;1;0;0;0;0;0;1 | |
1;0;0;0;0;0;0;1;0 | |
1;0;0;0;0;0;0;1;0 | |
1;0;0;0;0;0;0;1;0 | |
1;0;0;0;0;0;0;1;0 | |
1;0;0;0;0;0;0;1;0 | |
1;0;0;0;0;0;0;0;1 | |
1;0;0;0;0;0;0;0;1 | |
1;0;0;0;0;0;0;0;1 | |
1;0;0;0;0;0;0;0;1 | |
0;1;1;1;1;1;1;0;1 | |
0;1;1;1;1;1;1;0;1 | |
0;1;1;1;1;1;0;0;1 | |
0;1;1;1;1;0;0;1;0 | |
0;1;1;1;1;0;0;0;1 | |
0;1;1;1;1;0;0;0;1 | |
0;1;1;1;0;0;0;1;0 | |
0;1;1;1;0;0;0;0;1 | |
0;1;1;0;1;1;0;0;1 | |
0;1;1;0;0;0;0;1;0 | |
0;1;1;0;0;0;0;1;0 | |
0;1;1;0;0;0;0;1;0 | |
0;1;1;0;0;0;0;1;0 | |
0;1;1;0;0;0;0;1;0 | |
0;1;1;0;0;0;0;1;0 | |
0;1;1;0;0;0;0;1;0 | |
0;1;1;0;0;0;0;0;1 | |
0;1;1;0;0;0;0;0;1 | |
0;1;1;0;0;0;0;0;1 | |
0;1;1;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;1;0 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;1;0;0;0;0;0;0;1 | |
0;0;1;1;1;1;1;0;1 | |
0;0;1;1;1;1;1;0;1 | |
0;0;1;1;1;1;0;1;0 | |
0;0;1;1;1;1;0;0;1 | |
0;0;1;1;1;0;0;1;0 | |
0;0;1;1;1;0;0;1;0 | |
0;0;1;1;1;0;0;1;0 | |
0;0;1;1;1;0;0;1;0 | |
0;0;1;1;1;0;0;0;1 | |
0;0;1;1;1;0;0;0;1 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;1;0 | |
0;0;1;1;0;0;0;0;1 | |
0;0;1;1;0;0;0;0;1 | |
0;0;1;1;0;0;0;0;1 | |
0;0;1;1;0;0;0;0;1 | |
0;0;1;0;1;1;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;1;0 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;1;0;0;0;0;0;1 | |
0;0;0;1;1;1;1;1;0 | |
0;0;0;1;1;1;1;1;0 | |
0;0;0;1;1;1;1;1;0 | |
0;0;0;1;1;1;1;1;0 | |
0;0;0;1;1;1;1;1;0 | |
0;0;0;1;1;1;1;1;0 | |
0;0;0;1;1;1;1;0;1 | |
0;0;0;1;1;1;1;0;1 | |
0;0;0;1;1;1;0;1;0 | |
0;0;0;1;1;1;0;1;0 | |
0;0;0;1;1;1;0;1;0 | |
0;0;0;1;1;1;0;0;1 | |
0;0;0;1;1;1;0;0;1 | |
0;0;0;1;1;1;0;0;1 | |
0;0;0;1;1;1;0;0;1 | |
0;0;0;1;1;0;0;1;0 | |
0;0;0;1;1;0;0;1;0 | |
0;0;0;1;1;0;0;1;0 | |
0;0;0;1;1;0;0;1;0 | |
0;0;0;1;1;0;0;1;0 | |
0;0;0;1;1;0;0;1;0 | |
0;0;0;1;1;0;0;0;1 | |
0;0;0;1;1;0;0;0;1 | |
0;0;0;1;1;0;0;0;1 | |
0;0;0;1;1;0;0;0;1 | |
0;0;0;1;1;0;0;0;1 | |
0;0;0;1;0;1;1;0;1 | |
0;0;0;1;0;0;1;1;0 | |
0;0;0;1;0;0;1;0;1 | |
0;0;0;1;0;0;0;1;0 | |
0;0;0;1;0;0;0;1;0 | |
0;0;0;1;0;0;0;1;0 | |
0;0;0;1;0;0;0;1;0 | |
0;0;0;1;0;0;0;1;0 | |
0;0;0;1;0;0;0;1;0 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;1;0;0;0;0;1 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;1;0 | |
0;0;0;0;1;1;1;0;1 | |
0;0;0;0;1;1;1;0;1 | |
0;0;0;0;1;1;1;0;1 | |
0;0;0;0;1;1;1;0;1 | |
0;0;0;0;1;1;1;0;1 | |
0;0;0;0;1;1;1;0;1 | |
0;0;0;0;1;1;0;1;0 | |
0;0;0;0;1;1;0;1;0 | |
0;0;0;0;1;1;0;1;0 | |
0;0;0;0;1;1;0;0;1 | |
0;0;0;0;1;1;0;0;1 | |
0;0;0;0;1;1;0;0;1 | |
0;0;0;0;1;1;0;0;1 | |
0;0;0;0;1;1;0;0;1 | |
0;0;0;0;1;1;0;0;1 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;1;0 | |
0;0;0;0;1;0;0;0;1 | |
0;0;0;0;1;0;0;0;1 | |
0;0;0;0;1;0;0;0;1 | |
0;0;0;0;1;0;0;0;1 | |
0;0;0;0;1;0;0;0;1 | |
0;0;0;0;1;0;0;0;1 | |
0;0;0;0;1;0;0;0;1 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;1;0 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;1;0;1 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;1;0 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;1;0;0;1 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;1;0 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 | |
0;0;0;0;0;0;1;0;1 |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
#--------------------- Single-state | |
# read in data | |
dipper <- read.csv("dat/dipper.csv", sep = ";") | |
head(dipper) | |
dim(dipper) | |
# Create a sex variable with 1 = female and 2 = male | |
dipper$sex <- ifelse(dipper$males == 1, 2, 1) | |
head(dipper) | |
# Delete the redundant columns | |
dipper <- dipper[, -c(8,9)] | |
dipper <- as.matrix(dipper) | |
head(dipper) | |
nind <- nrow(dipper) | |
dipper.ch <- dipper[1:nind, 1:7] | |
dipper.sex <- dipper[1:nind, 8] | |
nocc <- ncol(dipper.ch) | |
# Create vector with occasion of marking | |
get.first <- function(x) min(which(x!=0)) | |
first <- apply(dipper.ch, 1, get.first) | |
# load nimble package | |
library(nimble) | |
# build model | |
code <- nimbleCode({ | |
# Priors | |
surv ~ dunif(0, 1) # Prior for survival | |
recap ~ dunif(0, 1) # Prior for recapture | |
for (i in 1:nind){ | |
# Define latent state at first capture | |
z[i,first[i]] <- 1 | |
for (t in (first[i]+1):nocc){ | |
# State process | |
z[i,t] ~ dbern(surv * z[i,t-1]) | |
# Observation process | |
y[i,t] ~ dbern(recap * z[i,t]) | |
} #t | |
} #i | |
}) | |
constants <- list(nind = nind, | |
nocc = nocc, | |
first = first) | |
data <- list(y = dipper.ch) | |
#--------- STANDARD APPROACH TO GENERATE INITS | |
# initial values for states | |
z.inits <- function(ch){ | |
state <- ch | |
state[state==0] <- 1 | |
get.first <- function(x) min(which(x!=0)) | |
f <- apply(ch, 1, get.first) | |
for (i in 1:nrow(ch)){ | |
state[i,1:f[i]] <- NA | |
} | |
return(state) | |
} | |
zInit <- z.inits(dipper.ch) | |
# initial values for all quantities | |
inits <- list(surv = 0.5, | |
recap = 0.9, | |
z = zInit) | |
# Create Model | |
Rmodel <- nimbleModel(code = code, | |
constants = constants, | |
data = data, | |
inits = inits, | |
check = FALSE, | |
calculate = FALSE) | |
Rmodel$surv | |
Rmodel$recap | |
Rmodel$z | |
Rmodel$calculate() | |
# Configure model | |
parameters <- c("surv", "recap", "z") | |
conf <- configureMCMC(Rmodel, | |
enableWAIC = TRUE, | |
monitors = parameters) | |
# Build model | |
Rmcmc <- buildMCMC(conf) | |
# Compile model | |
Cmodel <- compileNimble(Rmodel) | |
# Compile MCMC | |
Cmcmc <- compileNimble(Rmcmc, project = Rmodel) | |
# Run model | |
mcmc.out <- runMCMC(Cmcmc, | |
niter = 5000, | |
nburnin = 2500, | |
nchains = 2, | |
thin = 1, | |
WAIC = TRUE, | |
inits = inits, | |
samplesAsCodaMCMC = TRUE) | |
summary(mcmc.out$samples[,'surv']) | |
summary(mcmc.out$samples[,'recap']) | |
mcmc.out$WAIC | |
#--------- NIMBLE APPROACH TO GENERATE INITS | |
# initial values for "easy" quantities | |
inits <- list(surv = 0.5, recap = 0.9) | |
# Create Model | |
Rmodel <- nimbleModel(code = code, | |
constants = constants, | |
data = data, | |
inits = inits, | |
check = FALSE, | |
calculate = FALSE) | |
Rmodel$surv | |
Rmodel$recap | |
Rmodel$z | |
#Rmodel$getDependencies(c("surv", "recap"), determOnly = TRUE) | |
# initial values for "complex" quantities, latent states here | |
simNodes <- 'z' | |
simNodeScalar <- Rmodel$expandNodeNames(simNodes) | |
allNodes <- Rmodel$getNodeNames() | |
nodesSorted <- allNodes[allNodes %in% simNodeScalar] | |
set.seed(1) # This makes the simulations here reproducible | |
for(n in nodesSorted) { | |
Rmodel$simulate(n) | |
depNodes <- Rmodel$getDependencies(n) | |
Rmodel$calculate(depNodes) | |
} | |
Rmodel$z # z[i, first[i]] should be NA, it is not | |
Rmodel$calculate() # -Inf due to issue at line above, I guess | |
# change the state at first encounter to be NAs doesn't change anything | |
for (i in 1:nind){ | |
Rmodel$z[i,first[i]] <- NA | |
} | |
Rmodel$calculate('z') | |
Rmodel$calculate() # | |
# Configure model | |
parameters <- c("surv", "recap", "z") | |
conf <- configureMCMC(Rmodel, | |
enableWAIC = TRUE, | |
monitors = parameters) | |
# Build model | |
Rmcmc <- buildMCMC(conf) | |
# Compile model | |
Cmodel <- compileNimble(Rmodel) | |
# Compile MCMC | |
Cmcmc <- compileNimble(Rmcmc, project = Rmodel) | |
# Run model | |
mcmc.out <- runMCMC(Cmcmc, | |
niter = 5000, | |
nburnin = 2500, | |
thin = 1, | |
nchains = 2, | |
WAIC = TRUE, | |
inits = inits, | |
samplesAsCodaMCMC = TRUE) | |
summary(mcmc.out$samples[,'surv']) | |
summary(mcmc.out$samples[,'recap']) | |
mcmc.out$WAIC | |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment