Question: looks like the order of variables are inconsistent

[jiweiqi]

Question on: https://github.com/TuringLang/Turing.jl/edit/master/docs/_tutorials/10_bayesian-differential-equations.md

The order of variables in the model is different from the chain statistics. This will induce inconsistency for final plot of estimated trajectories.

Turing.setadbackend(:forwarddiff)

@model function fitlv(data, prob1)
    σ ~ InverseGamma(2, 3) # ~ is the tilde character
    α ~ truncated(Normal(1.5,0.5),0.5,2.5)
    β ~ truncated(Normal(1.2,0.5),0,2)
    γ ~ truncated(Normal(3.0,0.5),1,4)
    δ ~ truncated(Normal(1.0,0.5),0,2)

    p = [α,β,γ,δ]
    prob = remake(prob1, p=p)
    predicted = solve(prob,Tsit5(),saveat=0.1)

    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end

model = fitlv(odedata, prob1)

# This next command runs 3 independent chains without using multithreading.
chain = mapreduce(c -> sample(model, NUTS(.65),1000), chainscat, 1:3)

and

Chains MCMC chain (1000×17×3 Array{Float64,3}):

Iterations        = 1:1000
Thinning interval = 1
Chains            = 1, 2, 3
Samples per chain = 1000
parameters        = α, β, γ, δ, σ
internals         = acceptance_rate, hamiltonian_energy, hamiltonian_energy
_error, is_accept, log_density, lp, max_hamiltonian_energy_error, n_steps, 
nom_step_size, numerical_error, step_size, tree_depth

Summary Statistics
  parameters      mean       std   naive_se      mcse        ess      rhat
      Symbol   Float64   Float64    Float64   Float64    Float64   Float64

           α    1.2120    0.4926     0.0090    0.0911     6.0574   11.8201
           β    1.0672    0.1639     0.0030    0.0098   288.9165    1.0360
           γ    2.8014    0.2742     0.0050    0.0208    46.5485    1.0915
           δ    0.9538    0.0904     0.0017    0.0056   124.1565    1.0487
           σ    1.2712    0.6531     0.0119    0.1205     6.0766    9.8804

Quantiles
  parameters      2.5%     25.0%     50.0%     75.0%     97.5%
      Symbol   Float64   Float64   Float64   Float64   Float64

           α    0.5013    0.5243    1.5228    1.5737    1.6590
           β    0.6738    1.0262    1.0815    1.1337    1.4402
           γ    2.0953    2.6893    2.8377    2.9643    3.2437
           δ    0.8076    0.8993    0.9435    0.9901    1.1873
           σ    0.7482    0.7989    0.8394    2.1115    2.3534

and

chain_array = Array(chain)
for k in 1:300
    resol = solve(remake(prob1,p=chain_array[rand(1:1500), 1:4]),Tsit5(),saveat=0.1)
    plot!(resol, alpha=0.1, color = "#BBBBBB", legend = false)
end
# display(pl)
plot!(sol1, w=1, legend = false)

[jiweiqi]

@ChrisRackauckas I am not sure how to submit bugfix solutions. This is not a big deal but annoying to spend a while to figure out the problem.

[rikhuijzer]

It seems that disabling sorting is implemented at the MCMCChains.jl side, but not yet in Turing.jl. See https://github.com/TuringLang/MCMCChains.jl/pull/160 for details.

[devmotion]

The variables are not sorted in recent Turing versions by default: https://github.com/TuringLang/Turing.jl/pull/1627

Which Turing version did you use?

[jiweiqi]

The variables are not sorted in recent Turing versions by default: #1627

Which Turing version did you use?

Mine is v0.18.0. Any my full st:

Status `~/.julia/environments/v1.5/Project.toml`
  [d1b5debe] Arrhenius v0.1.1 `https://github.com/DENG-MIT/Arrhenius.jl#main`
  [c52e3926] Atom v0.12.30
  [fbb218c0] BSON v0.3.3
  [aae01518] BandedMatrices v0.16.10
  [052768ef] CUDA v2.4.3
  [39dd38d3] Dierckx v0.5.1
  [aae7a2af] DiffEqFlux v1.39.0
  [41bf760c] DiffEqSensitivity v6.49.1
  [0c46a032] DifferentialEquations v6.18.0
  [31c24e10] Distributions v0.25.14
  [587475ba] Flux v0.12.1
  [f6369f11] ForwardDiff v0.10.18
  [f67ccb44] HDF5 v0.15.6
  [a98d9a8b] Interpolations v0.13.2
  [42fd0dbc] IterativeSolvers v0.9.1
  [e5e0dc1b] Juno v0.8.4
  [a5e1c1ea] LatinHypercubeSampling v1.8.0
  [c7f686f2] MCMCChains v5.0.1
  [cc2ba9b6] MLDataUtils v0.5.4
  [429524aa] Optim v1.4.0
  [1dea7af3] OrdinaryDiffEq v5.60.1
  [91a5bcdd] Plots v1.19.4
  [49802e3a] ProgressBars v1.3.0
  [d330b81b] PyPlot v2.9.0
  [47a9eef4] SparseDiffTools v1.13.2
  [684fba80] SparsityDetection v0.3.4
  [f3b207a7] StatsPlots v0.14.28
  [c3572dad] Sundials v4.5.3
  [fce5fe82] Turing v0.18.0
  [fdbf4ff8] XLSX v0.7.6
  [ddb6d928] YAML v0.4.7
  [e88e6eb3] Zygote v0.6.12
  [8bb1440f] DelimitedFiles
  [37e2e46d] LinearAlgebra
  [de0858da] Printf
  [9a3f8284] Random
  [10745b16] Statistics

[devmotion]

Yeah but in the original comment above you just copied the output from the webpage (https://turing.ml/dev/tutorials/10-bayesian-differential-equations/), didn't you? And there a quite old version of Turing was used (0.15.18 apparently) that still sorted the variables. With recent versions of Turing the variables are not sorted anymore so you will get the order of variables in the chain that you expected.

I'm still not quite sure if and what the actual problem is here. Prior to the PR that I linked above (i.e. Turing < 0.16) the variables were sorted with natural sorting. So the order of the variables in the chain was expected and the desired behaviour. And hence also the indexing when plotting the trajectories seems to be correct.

[jiweiqi]

Yeah, everything I posted was from the document webpage. The problem is that I am using the new version of Turing.jl which will face the error of no-sorting. So the solution is we shall update the document according to the new version of Turing.jl.

[devmotion]

As explained in https://github.com/TuringLang/Turing.jl/pull/1627, you can recover the old behaviour (sorting with natural sort order) by specifying sort_chain=true (i.e., use sample(....; sort_chain=true)). Alternatively, you can just use 2:5 instead of 1:4 in the plotting code.

We should update the tutorial at some point but it is not completely urgent - in the tutorials the Turing version is stated on purpose so that you can see which version you have to use to get the same results. If you use a newer version, and in particular one with a breaking change, it is not guaranteed that the code still works.

[yebai]

Seems answered by https://github.com/TuringLang/Turing.jl/issues/1712#issuecomment-931624117. Please reopen if you have more questions.