Commit fad23301 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

solve problem in weight for subsampling and add init for stan

parent f9504b7e
### test stan
## I changed the coding of the random effect with norm(0,10*sigma it seems more slow!
## PROBLEM WITH STANDARDIZED VARIABLE DO I NEED TO INCLUDE AN INTERCEPT??
## VERY BAD CONVERGENCE ???
## Need to compare the two code
source('R/analysis/stan.run.R')
### TEST simple model on France only
df.lmer <- load.data.for.lmer(trait = 'SLA',fname = 'data.all.no.log.rds', cat.TF = FALSE,
sample.size = NA,
var.sample = NA,
sample.vec.TF. = FALSE,
select.set. = NA)
stan.list <- fun.turn.in.list.for.jags.stan(df.lmer,
cat.TF = FALSE)
source('R/analysis/model.stan/model.stan.LOGLIN.size.fixed.R', local = TRUE)
model <- load.model()
fun.init.stan <- function(chain_id= 1, stan.list){
init <- list(intercept = 0+(chain_id-1)/10,
intercept_species = rnorm(stan.list$N_species, 0, 0.5+(chain_id-1)/10),
intercept_plot = rnorm(stan.list$N_plot, 0, 0.5 +(chain_id-1)/10),
mean_logD = 2/3+(chain_id-1)/10,
sigma_inter_species = 0.5 +(chain_id-1)/10,
sigma_inter_plot = 0.5 + +(chain_id-1)/10,
sigma = 0.5 +(chain_id-1)/10)
}
inits <- list(model$init.fun(1, stan.list),
model$init.fun(2, stan.list),
model$init.fun(3, stan.list))
library(rstan)
system.time( stan.output <- stan(model_code = model$stan,
data = stan.list,
pars = model$pars,
init = inits,
iter = 1000,
warmup = 500,
chains = 3,
verbose = FALSE))
source('R/analysis/model.stan/model.stan.LOGLIN.ER.AD.Tf.fixed.oneset.R', local = TRUE)
model <- load.model()
fun.init.stan <- function(chain_id= 1, stan.list){
step <- (chain_id-1)/10
init <- list(intercept = 0+step,
intercept_species = rnorm(stan.list$N_species, 0, 0.5+step),
param_sumBn_species = rnorm(stan.list$N_species, 0, 0.5+step),
param_logD_species = rnorm(stan.list$N_species, 0, 0.5+step),
intercept_plot = rnorm(stan.list$N_plot, 0, 0.5 +step),
mean_logD = 2/3+step,
param_Tf = 2/3+step,
param_sumBn = 2/3+step,
param_sumTfBn = 2/3+step,
param_sumTnBn = 2/3+step,
param_sumTnTfBn_abs= 2/3+step,
sigma_inter_species = 0.5 +step,
sigma_inter_plot = 0.5 + +step,
sigma = 0.5 +step,
sigma_sumBn_species = 0.5 +step,
sigma_logD_species = 0.5 +step
)
}
inits <- list(fun.init.stan(1, stan.list),
fun.init.stan(2, stan.list),
fun.init.stan(3, stan.list))
library(rstan)
system.time( stan.output2 <- stan(model_code = model$stan,
data = stan.list,
pars = model$pars,
init = inits,
iter = 1000,
warmup = 500,
chains = 3,
verbose = FALSE))
save.image('test.stan.Rdata')
stan.output
plot(stan.output)
traceplot(stan.output)
pairs(stan.output)
library(ggmcmc)
S <- ggs(stan.output)
extract(stan.output, permuted = FALSE, inc_warmup=FALSE))
ggs_crosscorrelation(S, absolute_scale = FALSE)
ggs_density(S)
ggs_traceplot(S)
ggs_running(S)
ggs_autocorrelation(S)
ggs_Rhat(S)
ggs_geweke(S)
ggs_caterpillar(S)
kunstler@grp2233-Latitude-E6420.12167:1409576364
\ No newline at end of file
......@@ -19,6 +19,8 @@ model.files.jags.Tf.2b <-
model.files.jags.Tf.0 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.fixed.R")
model.files.jags.Tf.00 <-
c("R/analysis/model.jags/model.LOGLIN.fixed.R")
model.files.jags.Tf.3 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.biomes.R")
......
......@@ -39,9 +39,9 @@ model.files.lmer.Tf.2 <-
model.files.lmer.Tf.3 <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.biomes.species.R")
model.files.lmer.Tf.3b <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species.R")
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.biomes.species.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.fixed.species.R")
model.files.lmer.Tf.4 <-
c("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.species.R")
......@@ -125,7 +125,10 @@ output.dir <- function (type , model, trait, set, type.filling) {
## add sample equivalent per ecocode
add.sampling.prob.by.var.sample <- function(df, var.sample){
vec.prob <- 1/table(df[[var.sample]])
return(df)
df1 <- data.frame(ecocode = names(vec.prob),
prob.sample = as.numeric(vec.prob))
df2 <- merge(df, df1, by = 'ecocode')
return(df2)
}
......@@ -136,18 +139,23 @@ load.data.for.lmer <- function(trait, cat.TF,
sample.size. ,
var.sample. = 'ecocode',
select.biome. = NA,
select.set. = NA,
sample.vec.TF. = FALSE){
if (fname == 'data.all.no.log.rds') type <- 'no.log'
if (fname == 'data.all.log.rds') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- readRDS(file = file.path(base.dir,paste('data', trait, type,'rds',
if (fname == 'data.all.no.log.rds') type <- 'no.log'
if (fname == 'data.all.log.rds') type <- 'log'
if (cat.TF) cat <- 'cat'
if (!cat.TF) cat <- 'no.cat'
df <- readRDS(file = file.path(base.dir,paste('data', trait, type,'rds',
sep = '.')))
df$sp.name[is.na(df$sp.name)] <- 'missing.sp'
if(!is.na(select.biome.)) {
df <- df[df$biomes == select.biome., ]
}
if(!is.na(sample.size.)){
df$sp.name[is.na(df$sp.name)] <- 'missing.sp'
browser()
if(!is.na(select.biome.)) {
df <- df[df$biomes == select.biome., ]
}
if(!is.na(select.set.)) {
df <- df[df$set == select.set., ]
}
if(!is.na(sample.size.)){
if(!sample.vec.TF.){
if(sample.size. >nrow(df)){ sample.size. <- nrow(df)}
df <- add.sampling.prob.by.var.sample(df, var.sample = var.sample.)
......@@ -168,7 +176,7 @@ if(!is.na(sample.size.)){
type.filling = type.filling,
cat.TF = cat.TF)
print(paste('sub-sampled by',var.sample., 'by loading sample.vec'))
}
}
}else{
res <- fun.format.data.for.lmer(df, trait,
type.filling = type.filling,
......@@ -178,7 +186,7 @@ return( res)
}
load.and.save.data.for.lmer <- function(trait,
min.obs= 10,
min.obs= 10,
type.filling = species,
fname = 'data.all.no.log.rds',
base.dir = "output/processed"){
......
......@@ -10,9 +10,9 @@ load.model <- function(){
'sigma_inter_set', 'sigma_inter_plot',
'sigma_logD_species',# 'sigma_inter_tree',
'sigma_sumBn_species',
'sigma', 'vec_Tf_biomes',
'vec_sumBn_biomes', 'vec_sumTfBn_biomes',
'vec_sumTnBn_biomes', 'vec_sumTnTfBn_abs_biomes'),
'sigma', 'vec_Tf',
'vec_sumBn', 'vec_sumTfBn',
'vec_sumTnBn', 'vec_sumTnTfBn_abs'),
bug =
"#######################################################
######################## growth model with JAGS #######
......@@ -24,11 +24,11 @@ theo_g[i] <- inter[i] + DD[i] +TTf[i] + SUMBBN[i] +SUMTFBBN[i] + SUMTNBBN[i] + S
inter[i] <- intercept + intercept_species[species_id[i]] + intercept_plot[plot_id[i]] + intercept_set[set_id[i]] # + intercept_tree[tree_id[i]]*tree_01[i]
DD[i] <- (mean_logD + vec_logD_species[species_id[i]])*logD[i]
TTf[i] <- vec_Tf_biomes*Tf[i]
SUMBBN[i] <-(vec_sumBn_biomes + vec_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- (vec_sumTfBn_biomes)*sumTfBn[i]
SUMTNBBN[i] <- (vec_sumTnBn_biomes)*sumTnBn[i]
SUMABSBBN[i] <- (vec_sumTnTfBn_abs_biomes)*sumTnTfBn_abs[i]
TTf[i] <- vec_Tf*Tf[i]
SUMBBN[i] <-(vec_sumBn + vec_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- (vec_sumTfBn)*sumTfBn[i]
SUMTNBBN[i] <- (vec_sumTnBn)*sumTnBn[i]
SUMABSBBN[i] <- (vec_sumTnTfBn_abs)*sumTnTfBn_abs[i]
}
################################################
......@@ -47,11 +47,6 @@ for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_plot)
}
## tree effect
## for (j in 1:N_tree)
## {
## intercept_tree[j] ~ dnorm(0,tau_inter_tree)T(-7, 7)
## }
## set effect
for (j in 1:N_set)
......@@ -60,11 +55,11 @@ for (j in 1:N_set)
}
## biomes
vec_Tf_biomes ~ dnorm(0, tau0)T(-10, 10)
vec_sumBn_biomes ~ dnorm(0, tau0)T(-10, 10)
vec_sumTfBn_biomes ~ dnorm(0 ,tau0)T(-10, 10)
vec_sumTnBn_biomes ~ dnorm(0 ,tau0)T(-10, 10)
vec_sumTnTfBn_abs_biomes ~ dnorm(0 ,tau0)T(-10, 10)
vec_Tf ~ dnorm(0, tau0)T(-20, 20)
vec_sumBn ~ dnorm(0, tau0)T(-20, 20)
vec_sumTfBn ~ dnorm(0 ,tau0)T(-20, 20)
vec_sumTnBn ~ dnorm(0 ,tau0)T(-20, 20)
vec_sumTnTfBn_abs ~ dnorm(0 ,tau0)T(-20, 20)
###############################################
########### Non-hierarchical parameters ########
......@@ -72,8 +67,8 @@ for (j in 1:N_set)
tau0 <- 1.0E-4
#
intercept ~ dnorm(0,tau0)T(-5, 5)
mean_logD ~ dnorm(0,tau0)T(-5, 5)
intercept ~ dnorm(0,tau0)T(-15, 15)
mean_logD ~ dnorm(0,tau0)T(-15, 15)
# variance error
tau <- pow(sigma,-2)
......@@ -93,12 +88,11 @@ sigma_sumBn_species ~ dunif(0.0001,5)
#transformed variables
mean_Tf <- mean(vec_Tf_biomes)
mean_sumBn <- mean(vec_sumBn_biomes)
mean_sumTfBn <- mean(vec_sumTfBn_biomes)
mean_sumTnBn <- mean(vec_sumTnBn_biomes)
mean_sumTnTfBn_abs <- mean(vec_sumTnTfBn_abs_biomes)
mean_Tf <- mean(vec_Tf)
mean_sumBn <- mean(vec_sumBn)
mean_sumTfBn <- mean(vec_sumTfBn)
mean_sumTnBn <- mean(vec_sumTnBn)
mean_sumTnTfBn_abs <- mean(vec_sumTnTfBn_abs)
} # End of the jags model
")
......
load.model <- function(){
list(
name = "LOGLIN.fixed",
## parameters to save
pars = c('intercept' ,
'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_logD_species',# 'sigma_inter_tree',
'sigma_sumBn_species',
'sigma'),
bug =
"#######################################################
######################## growth model with JAGS #######
model {
############ Likelihood ###################
for (i in 1:N_indiv) {
logG[i] ~ dnorm(theo_g[i],tau)
theo_g[i] <- inter[i] + DD[i] + SUMBBN[i]
inter[i] <- intercept + intercept_species[species_id[i]] + intercept_plot[plot_id[i]] + intercept_set[set_id[i]]
DD[i] <- (mean_logD + vec_logD_species[species_id[i]])*logD[i]
SUMBBN[i] <-(vec_sumBn_species[species_id[i]])*sumBn[i]
}
################################################
########### Hierarchical parameters ########
### species random param
for (n in 1:N_species)
{
vec_logD_species[n] ~ dnorm(0,tau_logD_species)
vec_sumBn_species[n] ~ dnorm(0,tau_sumBn_species)
intercept_species[n] ~ dnorm(0,tau_inter_species)
}
## plot effect
for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_plot)
}
## set effect
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
tau0 <- 1.0E-4
#
intercept ~ dnorm(0,tau0)T(-15, 15)
mean_logD ~ dnorm(0,tau0)T(-15, 15)
# variance error
tau <- pow(sigma,-2)
sigma ~ dunif(0.00001,5)
tau_inter_plot <- pow(sigma_inter_plot,-2)
sigma_inter_plot ~ dunif(0.0001,5)
tau_inter_set <- pow(sigma_inter_set,-2)
sigma_inter_set ~ dunif(0.0001,5)
tau_inter_species <- pow(sigma_inter_species,-2)
sigma_inter_species ~ dunif(0.0001,5)
tau_logD_species <- pow(sigma_logD_species,-2)
sigma_logD_species ~ dunif(0.0001,5)
tau_sumBn_species <- pow(sigma_sumBn_species,-2)
sigma_sumBn_species ~ dunif(0.0001,5)
} # End of the jags model
")
}
load.model <- function(){
list(
name = "LOGLIN.fixed",
## parameters to save
pars = c('intercept' ,
'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_logD_species',# 'sigma_inter_tree',
'sigma_sumBn_species',
'sigma', 'param_sumBn'),
bug =
"#######################################################
######################## growth model with JAGS #######
model {
############ Likelihood ###################
for (i in 1:N_indiv) {
logG[i] ~ dnorm(theo_g[i],tau)
theo_g[i] <- inter[i] + DD[i] + SUMBBN[i]
inter[i] <- intercept + intercept_species[species_id[i]] + intercept_plot[plot_id[i]] + intercept_set[set_id[i]]
DD[i] <- (mean_logD + vec_logD_species[species_id[i]])*logD[i]
SUMBBN[i] <-(param_sumBn + vec_sumBn_species[species_id[i]])*sumBn[i]
}
################################################
########### Hierarchical parameters ########
### species random param
for (n in 1:N_species)
{
vec_logD_species[n] ~ dnorm(0,tau_logD_species)
vec_sumBn_species[n] ~ dnorm(0,tau_sumBn_species)
intercept_species[n] ~ dnorm(0,tau_inter_species)
}
## plot effect
for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_plot)
}
## set effect
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
tau0 <- 1.0E-4
#
intercept ~ dnorm(0,tau0)T(-15, 15)
mean_logD ~ dnorm(0,tau0)T(-15, 15)
#
param_sumBn ~ dnorm(0, tau0)T(-20, 20)
# variance error
tau <- pow(sigma,-2)
sigma ~ dunif(0.00001,5)
tau_inter_plot <- pow(sigma_inter_plot,-2)
sigma_inter_plot ~ dunif(0.0001,5)
tau_inter_set <- pow(sigma_inter_set,-2)
sigma_inter_set ~ dunif(0.0001,5)
tau_inter_species <- pow(sigma_inter_species,-2)
sigma_inter_species ~ dunif(0.0001,5)
tau_logD_species <- pow(sigma_logD_species,-2)
sigma_logD_species ~ dunif(0.0001,5)
tau_sumBn_species <- pow(sigma_sumBn_species,-2)
sigma_sumBn_species ~ dunif(0.0001,5)
} # End of the jags model
")
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.biomes.species",
list(name="lmer.LOGLIN.ER.AD.Tf.fixed.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+biomes.id+Tf+Tf:biomes.id+logD+sumBn+sumBn:biomes.id+sumTfBn+sumTfBn:biomes.id+sumTnBn+sumTnBn:biomes.id+sumTnTfBn.abs+sumTnTfBn.abs:biomes.id+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.fixed.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+biomes.id+Tf+Tf:biomes.id+logD+sumBn+sumBn:biomes.id+sumTfBn+sumTfBn:biomes.id+sumTnBn+sumTnBn:biomes.id+sumTnTfBn.abs+sumTnTfBn.abs:biomes.id+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)))
}
load.model <- function(){
list(name = "LOGLIN.ER.AD.Tf.r.set",
## parameters to save
list(name = "LOGLIN.ER.AD.Tf.fixed",
pars = c('intercept' , 'mean_logD', 'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_logD_species', 'sigma_sumBn_species', 'sigma',
'param_Tf', 'param_sumBn', 'param_sumTfBn',
'param_sumTnBn', 'param_sumTnTfBn_abs'),
## parameters to save
init.fun = function(chain_id= 1, list){
step <- (chain_id-1)/10
init <-
list(intercept = 0+step,
intercept_species = rnorm(list$N_species, 0, 0.5+step),
param_sumBn_species = rnorm(list$N_species, 0, 0.5+step),
param_logD_species = rnorm(list$N_species, 0, 0.5+step),
intercept_plot = rnorm(list$N_plot, 0, 0.5 +step),
intercept_set = rnorm(list$N_set, 0, 0.5 +step),
mean_logD = 2/3+step,
param_Tf = 2/3+step,
param_sumBn = 2/3+step,
param_sumTfBn = 2/3+step,
param_sumTnBn = 2/3+step,
param_sumTnTfBn_abs= 2/3+step,
sigma_inter_species = 0.5 +step,
sigma_inter_plot = 0.5 + +step,
sigma_inter_set = 0.5 + +step,
sigma = 0.5 +step,
sigma_sumBn_species = 0.5 +step,
sigma_logD_species = 0.5 +step)
return(init)
},
stan = "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_set;
int<lower=0> N_tree;
int species_id[N_indiv];
int plot_id[N_indiv];
int set_id[N_indiv];
int tree_id[N_indiv];
real logG[N_indiv];
real logD[N_indiv];
real sumBn[N_indiv];
......@@ -24,40 +51,27 @@ data {
real sumTfBn[N_indiv];
real sumTnBn[N_indiv];
real sumTnTfBn_abs[N_indiv];
real tree_01[N_indiv];
}
parameters{
real<lower=-1,upper=1> intercept;
real<lower=-10,upper=10> mean_logD;
real<lower=-10,upper=10> mean_Tf;
real<lower=-10,upper=10> mean_sumBn;
real<lower=-10,upper=10> mean_sumTfBn;
real<lower=-10,upper=10> mean_sumTnBn;
real<lower=-10,upper=10> mean_sumTnTfBn_abs;
real<lower=-5,upper=5> intercept;
real<lower=-5,upper=5> mean_logD;
real<lower=0.0001,upper=5> sigma_inter_species;
real<lower=0.0001,upper=5> sigma_inter_set;
real<lower=0.0001,upper=5> sigma_inter_plot;
real<lower=0.0001,upper=5> sigma_inter_tree;
real<lower=0.0001,upper=5> sigma_logD_species;
real<lower=0.0001,upper=5> sigma_Tf_set;
real<lower=0.0001,upper=5> sigma_sumBn_set;
real<lower=0.0001,upper=5> sigma_sumBn_species;
real<lower=0.0001,upper=5> sigma_sumTfBn_set;
real<lower=0.0001,upper=5> sigma_sumTnBn_set;
real<lower=0.0001,upper=5> sigma_sumTnTfBn_abs_set;
real<lower=0.0001,upper=5> sigma;
vector[N_set] param_Tf_set;
vector[N_set] param_sumBn_set;
vector[N_species] param_sumBn_species;
vector[N_set] param_sumTfBn_set;
vector[N_set] param_sumTnBn_set;
vector[N_set] param_sumTnTfBn_abs_set;
vector[N_species] param_logD_species;
vector[N_species] intercept_species;
vector[N_plot] intercept_plot;
vector[N_set] intercept_set;
vector[N_tree] intercept_tree;
real<lower=-15,upper=15> param_Tf;
real<lower=-15,upper=15> param_sumBn;
real<lower=-15,upper=15> param_sumBn_species[N_species];
real<lower=-15,upper=15> param_sumTfBn;
real<lower=-15,upper=15> param_sumTnBn;
real<lower=-15,upper=15> param_sumTnTfBn_abs;
real<lower=-15,upper=15> param_logD_species[N_species];
real<lower=-15,upper=15> intercept_species[N_species];
real<lower=-15,upper=15> intercept_plot[N_plot];
real<lower=-15,upper=15> intercept_set[N_set];
}
transformed parameters {
# vector for prediction
......@@ -65,22 +79,21 @@ transformed parameters {
for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_species[species_id[i]]
+ intercept_tree[tree_id[i]]*tree_01[i]
+ intercept_plot[plot_id[i]]
+ intercept_set[set_id[i]]
+ (mean_logD + param_logD_species[species_id[i]])*logD[i]
+ (mean_Tf + param_Tf_set[set_id[i]])*Tf[i]
+ (mean_sumBn + param_sumBn_set[set_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
+ (mean_sumTfBn + param_sumTfBn_set[set_id[i]])*sumTfBn[i]
+ (mean_sumTnBn + param_sumTnBn_set[set_id[i]])*sumTnBn[i]
+ (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_set[set_id[i]])*sumTnTfBn_abs[i]
;
+ (mean_logD + param_logD_species[species_id[i]])*logD[i]
+ (param_Tf)*Tf[i]
+ ( param_sumBn + param_sumBn_species[species_id[i]])*sumBn[i]
+ (param_sumTfBn)*sumTfBn[i]
+ (param_sumTnBn)*sumTnBn[i]
+ (param_sumTnTfBn_abs)*sumTnTfBn_abs[i]
;
}
}
model {
# constants for prior