Skip to content

Instantly share code, notes, and snippets.

@oliviergimenez

oliviergimenez/dipper.csv

Created December 3, 2020 04:41
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
  • Save oliviergimenez/e41e9cb99174f2124f948308e19ca7ec to your computer and use it in GitHub Desktop.
Save oliviergimenez/e41e9cb99174f2124f948308e19ca7ec to your computer and use it in GitHub Desktop.
Download ZIP
Use Nimble to generate inits for latent states
Raw
dipper.csv
We can make this file beautiful and searchable if this error is corrected: No commas found in this CSV file in line 0.
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
Raw
generate-inits-ssm.R
#--------------------- 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