Chains are not independent

[itsdfish]

Hi Rob-

I noticed that the chains are not independent. Another potential issue is that the functionality of ncommands and nchains is unclear. Both appear to run n instances of the same chain, but ncommands appears to run them in parallel. Is this a bug or a problem with my setup?

Here is a simple example:

ENV["JAGS_HOME"] = "usr/bin/jags" #your path here
using Jags, StatsPlots, Random, Distributions
#cd(@__DIR__)
ProjDir = pwd()
Random.seed!(3431)

y = rand(Normal(0,1),50)

Model = "
model {
      for (i in 1:length(y)) {
            y[i] ~ dnorm(mu,sigma);
      }
      mu  ~ dnorm(0, 1/sqrt(10));
      sigma  ~ dt(0,1,1) T(0, );
  }
"

monitors = Dict(
  "mu" => true,
  "sigma" => true,
  )

jagsmodel = Jagsmodel(
  name="Gaussian",
  model=Model ,
  monitor=monitors,
  ncommands=4, nchains=1,
  #deviance=true, dic=true, popt=true,
  pdir=ProjDir
  )

println("\nJagsmodel that will be used:")
jagsmodel |> display

data = Dict{String, Any}(
  "y" => y,
)

inits = [
  Dict("mu" => 0.0,"sigma" => 1.0,
  ".RNG.name" => "base::Mersenne-Twister")
]

println("Input observed data dictionary:")
data |> display
println("\nInput initial values dictionary:")
inits |> display
println()
#######################################################################################
#                                 Estimate Parameters
#######################################################################################
sim = jags(jagsmodel, data, inits, ProjDir)
sim = sim[5001:end,:,:]
plot(sim)

[goedman]

Hi Chris, why do you think that? I see 4 chains in the plots. In the line1 example I don't set the seed.

[goedman]

nchains is the total number of chains created, ncommands set the number of processes used to run the nchains. Not sure how I handled the case where commands > 1 and not equal to chains.

[itsdfish]

When I set nchains = ncommands = 4, it provides 16 identical chains.

sim.value[:,:,1] == sim.value[:,:,2]
true

For visual confirmation, the plot below shows that all chains are superimposed. Are you not getting this result on your machine?

[itsdfish]

Maybe I introduced a bug in one of my updates

[goedman]

Hmm,

Jagsmodel that will be used:
  name =                    "line1"
  ncommands =               4
  nchains =                 2
  adapt =                   1000
  nsamples =                10000
  thin =                    1
  jagsthin =                1
  monitor =                 Dict{String,Bool}("sigma" => 1,"alpha" => 1,"tau" => 1,"beta" => 1)
  deviance =                true
  dic =                     true
  popt =                    true
  model_file =              line1.bugs
  data_file =               line1-data.R
  tmpdir =                  /Users/rob/.julia/dev/Jags/Examples/Line1/tmp

Input observed data dictionary:
Dict{String,Any} with 3 entries:
  "n" => 5
  "x" => [1, 2, 3, 4, 5]
  "y" => [1, 3, 3, 3, 5]

Input initial values dictionary:
1-element Array{Dict{String,Int64},1}:
 Dict("alpha" => 0,"tau" => 1,"beta" => 0)


Creating data file line1-data.R -   0.000638 seconds (100 allocations: 5.375 KiB)

Creating init file line1-inits1.R -   0.000562 seconds (44 allocations: 2.891 KiB)
Creating init file line1-inits2.R -   0.000465 seconds (44 allocations: 2.891 KiB)
Creating init file line1-inits3.R -   0.000507 seconds (44 allocations: 2.891 KiB)
Creating init file line1-inits4.R -   0.000597 seconds (44 allocations: 2.891 KiB)

Executing 4 command(s), each with 2 chain(s) took:
  0.288963 seconds (133 allocations: 5.703 KiB)

Reading line1-cmd1-chain1.txt
Reading line1-cmd1-chain2.txt
Reading line1-cmd2-chain1.txt
Reading line1-cmd2-chain2.txt
Reading line1-cmd3-chain1.txt
Reading line1-cmd3-chain2.txt
Reading line1-cmd4-chain1.txt
Reading line1-cmd4-chain2.txt

Object of type Chains, with data of type 10000×5×8 Array{Float64,3}

Iterations        = 1:10000
Thinning interval = 1
Chains            = 1, 2, 3, 4, 5, 6, 7, 8
Samples per chain = 10000
parameters        = alpha, beta, deviance, sigma, tau

2-element Array{ChainDataFrame,1}

Summary Statistics

│ Row │ parameters │ mean     │ std      │ naive_se   │ mcse       │ ess     │ r_hat    │
│     │ Symbol     │ Float64  │ Float64  │ Float64    │ Float64    │ Any     │ Any      │
├─────┼────────────┼──────────┼──────────┼────────────┼────────────┼─────────┼──────────┤
│ 1   │ alpha      │ 3.00139  │ 0.574499 │ 0.00203116 │ 0.00232841 │ 78327.8 │ 1.00011  │
│ 2   │ beta       │ 0.799041 │ 0.408674 │ 0.00144488 │ 0.00153812 │ 69870.1 │ 1.00028  │
│ 3   │ deviance   │ 12.9037  │ 3.69719  │ 0.0130715  │ 0.0252219  │ 20302.5 │ 0.999944 │
│ 4   │ sigma      │ 1.01157  │ 0.761745 │ 0.00269318 │ 0.00590691 │ 16194.8 │ 0.999991 │
│ 5   │ tau        │ 1.87919  │ 1.53463  │ 0.00542575 │ 0.00798475 │ 39459.0 │ 0.999944 │

Quantiles

│ Row │ parameters │ 2.5%     │ 25.0%    │ 50.0%    │ 75.0%    │ 97.5%   │
│     │ Symbol     │ Float64  │ Float64  │ Float64  │ Float64  │ Float64 │
├─────┼────────────┼──────────┼──────────┼──────────┼──────────┼─────────┤
│ 1   │ alpha      │ 1.96614  │ 2.75586  │ 3.00052  │ 3.25033  │ 4.02157 │
│ 2   │ beta       │ 0.052762 │ 0.621173 │ 0.800534 │ 0.978123 │ 1.52614 │
│ 3   │ deviance   │ 8.79272  │ 10.241   │ 11.9244  │ 14.5078  │ 22.343  │
│ 4   │ sigma      │ 0.412698 │ 0.625137 │ 0.817692 │ 1.14762  │ 2.71555 │
│ 5   │ tau        │ 0.135607 │ 0.759282 │ 1.49562  │ 2.55888  │ 5.87132 │


julia> sim[:,:,1] == sim[:,:,5]
false

julia> sim[:,1,1] == sim[:,1,2]
false

julia> sim[:,1,1] == sim[:,1,3]
false

julia> sim[:,1,1] == sim[:,1,4]
false

julia> sim[:,1,1] == sim[:,1,5]
false

[goedman]

Just ncommands=1, nchains=2:

GKSTerm.pdf

I might have to check if tmp is cleared properly. My statement above is wrong, 4 chains in parallel is ncommands=4, nchains=1.

[itsdfish]

Strange. Do you also obtain the correct results with the code in my initial post?

[itsdfish]

P.S. I think you need to test sim.value[:,:,1] == sim.value[:,:,2] instead of the chain. On my system I get:

julia> sim.value[:,:,1] == sim.value[:,:,2]
true

julia> sim[:,:,1] == sim[:,:,2]
false

julia> sim.value[1:5,:,1]
2-dimensional AxisArray{Float64,2,...} with axes:
    :iter, 5001:1:5005
    :var, ["mu", "sigma"]
And data, a 5×2 Array{Float64,2}:
  0.215479   0.705492
  0.26205    0.6893  
  0.0842925  0.792049
 -0.0962303  0.662537
 -0.0810468  0.805239

julia> sim.value[1:5,:,2]
2-dimensional AxisArray{Float64,2,...} with axes:
    :iter, 5001:1:5005
    :var, ["mu", "sigma"]
And data, a 5×2 Array{Float64,2}:
  0.215479   0.705492
  0.26205    0.6893  
  0.0842925  0.792049
 -0.0962303  0.662537
 -0.0810468  0.805239

[itsdfish]

I tried to compare version 2.0.1 to the most current version but I couldn't get it to work any more. Really strange.

In case I set up my example incorrectly, could you try it and/or send me the example that was successful for you? I'll give it a try once I get back from the gym.

[itsdfish]

Ok. After finding your example and comparing it to mine, I discovered part of the problem. The code I copied and pasted to create my model had a random seed in the inits dictionary:

inits = [
  Dict("mu" => 0.0,"sigma" => 1.0,
  ".RNG.name" => "base::Mersenne-Twister")
]

After removing the RNG variable, I was able to get different chains in some cases, but not others.
With the setup below, I get 4 identical chains in parallel:

jagsmodel = Jagsmodel(
  name="Gaussian",
  model=Model ,
  monitor=monitors,
  ncommands=4, nchains=1,
  update = 100000,
  #deviance=true, dic=true, popt=true,
  pdir=ProjDir
  )

With the setup below, the result is 4 different chains sampled sequentially:

jagsmodel = Jagsmodel(
  name="Gaussian",
  model=Model ,
  monitor=monitors,
  ncommands=1, nchains=4,
  update = 100000,
  #deviance=true, dic=true, popt=true,
  pdir=ProjDir
  )

Finally, the setup below yields 4 different sets of chains and chains with each set are identical.

jagsmodel = Jagsmodel(
  name="Gaussian",
  model=Model ,
  monitor=monitors,
  ncommands=4, nchains=4,
  update = 100000,
  #deviance=true, dic=true, popt=true,
  pdir=ProjDir
  )

Confirming this is the case:

julia> sum(map(x->sim.value[:,:,1] == sim.value[:,:,x],1:16))
4

So is there a way to run 4 different chains in parallel?

[goedman]

Not right now, looking at the code, the way to achieve this is to pass a cmd-x-chain-y specific seed in the Gaussian-cmdx.jags files, e.g.:

/*
	Generated Gaussian.jags command file
*/
model in Gaussian.bugs
data in Gaussian-data.R
compile, nchains(2)
parameters in Gaussian-inits2.R, chain(1)
parameters in Gaussian-inits1.R, chain(2)
initialize
update 1000
monitor sigma, thin(1)
monitor mu, thin(1)
update 2000
coda *, stem(Gaussian-cmd2-)
exit

needs to be expanded to:

/*
	Generated Gaussian.jags command file
*/
model in Gaussian.bugs
data in Gaussian-data.R
compile, nchains(2)
parameters in Gaussian-cmd1-inits2.R, chain(1)
parameters in Gaussian-cmd1-inits1.R, chain(2)
initialize
update 1000
monitor sigma, thin(1)
monitor mu, thin(1)
update 2000
coda *, stem(Gaussian-cmd2-)
exit

for cmd1. Etc.This might be a consequence of the fact this needs to be compiled in.

Definitely a bit of work. I don't think that has ever worked properly in earlier releases. Shows that Jags.jl hasn't been used that much. How important is it for your work?

[goedman]

To documents this a bit further, each cmd and each chain need an init file generated from:

inits = [
  Dict("mu" => 0.0,"sigma" => 1.0,
    ".RNG.name" => "base::Mersenne-Twister",
    ".RNG.seed" => the-cmd-x-chain-y-specific-seed
  )
]

in your example.

[goedman]

And the jagsmodel.command array:

julia> jagsmodel.command
2-element Array{Base.AbstractCmd,1}:
 `jags Gaussian-cmd1.jags`
 `jags Gaussian-cmd2.jags`

in the case of 4 parallel executions, each with 2 chains needs to look like:

julia> jagsmodel.command
8-element Array{Base.AbstractCmd,1}:
 `jags Gaussian-cmd1-chain1.jags`
 `jags Gaussian-cmd1-chain2.jags`
 `jags Gaussian-cmd2-chain1.jags`
 `jags Gaussian-cmd2-chain2.jags`
...
 `jags Gaussian-cmd4-chain2.jags`

[itsdfish]

Thanks, Rob. Its somewhat important. I started to use Jags because Turing does not currently have an efficient sampler for discrete parameters (see this issue). They may have a solution in 3-6 months, but it's often the case that projects take longer than anticipated.

Do you think adding this capability would be difficult?

[goedman]

Not that difficult, but just a very detail-sensitive project. Basically something like:

1. Generate ncommands x nchains seeds (this needs to be reproducible).
2. Generate input files for all init files (including entries for these seed values) in the tmp dir
3. Generate appropriate name-cmd-x-chain-y.jags files
4. Create the model.command array

I think those are the steps needed.

[itsdfish]

Hi Rob,

I just wanted to follow up with you about this issue. As a summary, Jags.jl runs identical chains in parallel, which defeats the purpose of having multiple parallel chains. Do you still have plans to fix this?

[goedman]

Hi Chris,

In the Jags docs, on page 16:

A suitable set of initial values (for a single chain) would be
list(alpha = 0, beta = c(NA, 1.03, -2.03, 0.52))
This allows parameter values to be supplied for the stochastic elements of beta but not the constant first element.
If initial values are not supplied by the user, then each parameter chooses its own initial value based on the values of its parents. The initial value is chosen to be a “typical value” from the prior distribution. The exact meaning of “typical value” depends on the distribution of the stochastic node, but is usually the mean, median, or mode.
If you rely on automatic initial value generation and are running multiple parallel chains, then the initial values will be the same in all chains.1 You are advised to set the starting values manually.

with a footnote:

This is undesirable behaviour and it will be changed in a future release of JAGS.

My conclusion is that:

...
jagsmodel = Jagsmodel(
  name="line1",
  model=line,
  ncommands=4,
  nchains=4,
  monitor=monitors,
  deviance=true, dic=true, popt=true,
  pdir=ProjDir
  );

println("\nJagsmodel that will be used:")
jagsmodel |> display

data = Dict(
  "x" => [1, 2, 3, 4, 5],
  "y" => [1, 3, 3, 3, 5],
  "n" => 5
)

inits = [
  Dict("alpha" => 0,"beta" => 0,"tau" => 1),
  Dict("alpha" => 1,"beta" => 2,"tau" => 1),
  Dict("alpha" => 3,"beta" => 3,"tau" => 2),
  Dict("alpha" => 5,"beta" => 2,"tau" => 5)
]
...

might achieve this.

But I do have a problem with the Jags binary in that ncommands in above script used to be 1, but now it hangs Julia if I don't set it to 4. This could also be the new thread model in Julia. I do see 4 chains.

[goedman]

If I recall correctly, my intention was that using chains Jags controls running the individual chains, with ncommands Julia is supposed to be the master.

[itsdfish]

Hi Rob,

Thank you for clarifying. Would it be possible to add your example above to the README as an example of running multiple independent chains? I think that will be useful for other people who use Jags.jl. Aside from that, I think the issue is resolved. Thanks again!

[goedman]

Hi Chris,

Yes, as far as possible, it is, but I'm worried about why Julia hangs with 4 Julia commands each with 1 chain.

Will try to find some time to dig into this. Edit: Still working on the revamp of StanJulia. There is now a StanSample.jl v5.0.0 branch which is all keyword argument based but before releasing it I need to do the others (StanOptimize and StanVariational).

[itsdfish]

Sorry. I overlooked that. Thank for looking into this issue.

Great news about StanSample.jl by the way. I will update as soon as it is released!

[goedman]

It's still a branch only (v5).

I actually wonder if one of the random number generators in Jags is failing (on MacOS) and this causes Julia to hang.

[goedman]

Doesn't seem to be related to setting different random number generators but it is the Jags binary hanging on my system.

Did put a note on multiple inits in the readme and a warning about the problem.

For now I will leave this issue open, not sure if I can fix it.

By the way, have you tried Mamba for Gibbs?

[itsdfish]

Thanks for looking into the issue. I played around with Mamba several years ago but did not think about using it for Gibbs. Thanks for hte idea.