Test code

This is some test code I wrote, that I want to keep in case I need to come back to it:

model_to_data_block("model { for(i in 1:n) {b[i] ~ 1}}") %>% cat

remove_priors("data{alpha ~ dunif(0, 1)\nbeta <- alpha}\n", c(alpha = 0.5)) %>% cat

set.seed(101) set_seed(list())

set.seed(102) generate_data("model{beta ~ dunif(0,up) \n for(i in 1:5){alpha[i,2] = 2}}", monitor = c("beta", "alpha"), inits = list(), parameters = list("up"=1), data = list())

model <- "model{

x[1:2] ~ dmnorm(mu[1:2], Omega[1:2,1:2])

mu[1] <- 0

mu[2] <- 0

Omega[1,1] <- 1

Omega[1,2] <- 0.9

Omega[2,1] <- 0.9

Omega[2,2] <- 1

a <- x[1]

b <- x[2]

}

"

#Code to simulate from a bivariate normal model <- "model{ X ~ dnorm(5,1) e ~ dnorm(0.5,1) Y <- 0.5*X + e } "

data <- list() monitor <- c("X", "Y") inits = list(list())

library(R2jags) set.seed(101) samples <- jags(data=data, inits=inits, parameters.to.save=monitor, model.file = textConnection(model), n.chains = 1, n.iter = 50000, n.burnin = 0, n.thin = 1, DIC=FALSE)

plot(samples$BUGSoutput$sims.list$X, samples$BUGSoutput$sims.list$Y) #looks like samples from a bivariate normal

acf(samples$BUGSoutput$sims.list$X, lag.max =5) #autocorrelation plot show no autocorrelation for X acf(samples$BUGSoutput$sims.list$Y, lag.max =5) #autocorrelation plot show no autocorrelation for Y

samples <- jags(data=list("e"=0), inits=inits, parameters.to.save=monitor, model.file = textConnection(model), n.chains = 1, n.iter = 50000, n.burnin = 0, n.thin = 1, DIC=FALSE)

plot(samples$BUGSoutput$sims.list$X, samples$BUGSoutput$sims.list$Y) #looks like samples from a bivariate normal

acf(samples$BUGSoutput$sims.list$X, lag.max =5) #autocorrelation plot show no autocorrelation for X acf(samples$BUGSoutput$sims.list$Y, lag.max =5) #autocorrelation plot show no autocorrelation for Y

#runjags samples2 <- run.jags(data=list("dummy"=0), monitor=monitor, model=model, n.chains=1, sample=50000, burnin=0, thin=1, summarise=FALSE) samples2$modules samples2$factories samples2 <- run.jags(data=list("X"=0), monitor=monitor, model=model, n.chains=1, sample=50000, burnin=0, thin=1, summarise=FALSE) samples2$modules samples2$factories

#https://appsilon.com/how-to-sample-from-multidimensional-distributions-using-gibbs-sampling/

gibbs_normal_sampling <- function(n_iteration, init_point, ro) {

init point is some numeric vector of length equals to 2

theta_1 <- c(init_point[1], numeric(n_iteration)) theta_2 <- c(init_point[2], numeric(n_iteration)) for (i in 2:(n_iteration+1)) { theta_1[i] <- rnorm(1, ro * theta_2[i-1], sqrt(1 - ro^2)) theta_2[i] <- rnorm(1, ro * theta_1[i], sqrt(1 - ro^2)) } list(theta_1, theta_2) }

te <- gibbs_normal_sampling(50000, c(0,0), 0.9) acf(te1)

library(runjags)

Likelihood: Note key components of a GLM on one line each

model <- "data{for (i in 1:n){ C[i] ~ dpois(lambda[i]) # 1. Distribution for random part log(lambda[i]) <- log.lambda[i] # 2. Link function log.lambda[i] <- alpha + beta1 * year[i] + beta2 * pow(year[i],2) + beta3 * pow(year[i],3) # 3. Linear predictor } #i bla = C[1] + C[2] } model{ dummy <- 0 }"

library(runjags) model <- "model{ X ~ dnorm(0,1) e ~ dnorm(0,1) Y <- 0.5*X + e } " samples2 <- run.jags(data = list(list()), monitor=c("X", "Y"), model = model, n.chains = 1, sample = 2, burnin = 0, thin = 1, summarise=FALSE)

names(samples2) samples2$end.state

library(rjags) rjags::list.factories("sampler") rjags::unload.module("bugs") rjags::list.factories("sampler") rjags::set.factory("base::Finite", "sampler", FALSE) rjags::set.factory("base::Slice", "sampler", FALSE) rjags::list.factories("sampler")

rjags::set.factory("base::Slice", "sampler", TRUE) rjags::list.modules() rjags::unload.module("glm") rjags::unload.module("bugs") rjags::list.modules()