Search code examples
bayesianjagsrjags

How to randomly initialize the chains with Rjags?


I'm trying initialize the chains of the bayesian model randomly, today i do it inserting the values manually, like this:

inits1 <- list("alpha" = 10, "beta" = 3.2, "lambda" = 2)
inits2 <- list("alpha" = 18, "beta" = 7, "lambda" = 4.6)
inits3 <- list("alpha" = 16, "beta" = 5, "lambda" = 3.1)

model.inits <- list(inits1, inits2, inits3)

I tried performed it by the the following way, but resulted in error:

inits <- function(){ list("alpha"=rnorm(p), "beta" = rlnorm(p), "lambda" = runif(p))}

Where p->3 ( the number of parameters of the model)

Follows the error:

Error: The following error occured when compiling and adapting the model using rjags: Error in setParameters(init.values[[i]], i) : RUNTIME ERROR: Dimension mismatch in values supplied for alpha


Solution

  • Here is an example of an initialization function for rjags that would work for up to a total of 8 chains:

    inits_func <- function(chain){
      gen_list <- function(chain = chain){
        list( 
          alpha = rnorm(1, 14, 2),
          beta = rlnorm(1, log(5), log(1)),
          lambda = runif(1, 0.05, 5),
          .RNG.name = switch(chain,
                             "1" = "base::Wichmann-Hill",
                             "2" = "base::Marsaglia-Multicarry",
                             "3" = "base::Super-Duper",
                             "4" = "base::Mersenne-Twister",
                             "5" = "base::Wichmann-Hill",
                             "6" = "base::Marsaglia-Multicarry",
                             "7" = "base::Super-Duper",
                             "8" = "base::Mersenne-Twister"),
          .RNG.seed = sample(1:1e+06, 1)
        )
      }
      return(switch(chain,           
                    "1" = gen_list(chain),
                    "2" = gen_list(chain),
                    "3" = gen_list(chain),
                    "4" = gen_list(chain),
                    "5" = gen_list(chain),
                    "6" = gen_list(chain),
                    "7" = gen_list(chain),
                    "8" = gen_list(chain)
      )
      )
    }
    

    The important thing to note is that the function should 1) have a single argument named chain and 2) you need to manually specify the number of parameters to get estimated for each chain.

    In your function you have specified p=3. Therefore you end up generating three alpha, beta, and lambda parameters per chain. Looking at the lists you used before, it would appear that each chain should have a single alpha, beta, and lambda parameters (i.e., there is only one alpha, beta, and lambda parameter in your model). Additionally, you should specify the parameters for (e.g., mu and sd for the normal distribution, the upper and lower bound for the uniform distribution). I 'eyeballed' these based off the values you specified in the three lists. You could greatly simplify this function as well and not specify the RNG algorithm's for each chain if you wanted to. The function you wrote could work if instead it was written as:

    inits <- function(chain){
     list(
     alpha = rnorm(1),
     beta = rlnorm(1),
     lambda = runif(1))
    }
    

    But again, you'd want to specify the other parameters in rnorm, rlnorm, and runif.