Order of the multi-level parameter is not preserved when converted into InferenceData

[AnselmJeong]

Thank you for this remarkable package...
As a newbie, I'm trying to replicate the centered_eight example from arviz homepage.

I've obtained successfully the posterior sample from the model. (let's call this chain chn_post)
And the I've created an InferenceData from this posterior by from_mcmcchains() function. (=idata)
This idata must contain the above mentioned chn_post as idata.posterior.

Below is the result of summarystats of the two version of the chain.

    summarystats(chn_post)

     parameters      mean       std   naive_se      mcse         ess      rhat    ⋯
   
           μ    4.3999    3.2838     0.0367    0.0976   1334.4223    1.0020    ⋯
           τ    4.3683    3.0589     0.0342    0.1179    713.4324    1.0097    ⋯
        θ[1]    6.6081    5.9446     0.0665    0.1327   2420.8891    1.0019    ⋯
        θ[2]    5.0276    4.9369     0.0552    0.0980   2592.4887    1.0014    ⋯
        θ[3]    3.7422    5.4718     0.0612    0.1177   2451.2714    1.0005    ⋯
        θ[4]    4.8221    5.0381     0.0563    0.1072   2825.2853    1.0010    ⋯
        θ[5]    3.4301    4.8281     0.0540    0.1178   2073.0898    1.0016    ⋯
        θ[6]    3.9433    5.0540     0.0565    0.1108   2378.3583    1.0015    ⋯
        θ[7]    6.7735    5.4439     0.0609    0.1182   2091.5826    1.0012    ⋯
        θ[8]    4.9105    5.6013     0.0626    0.1054   3117.7890    1.0010    ⋯

    summarystats(idata.posterior)

     variable	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail

	μ	4.4	3.284	-1.781	10.603	0.09	0.064	1336.0	1866.0
	τ	4.368	3.059	0.936	9.741	0.114	0.081	366.0	175.0
	θ[1]	3.43	4.828	-6.15	12.045	0.106	0.075	1995.0	3179.0
	θ[2]	4.822	5.038	-5.274	13.899	0.095	0.067	2543.0	2929.0
	θ[3]	4.91	5.601	-5.469	15.622	0.1	0.074	2676.0	3293.0
	θ[4]	3.943	5.054	-6.119	13.205	0.104	0.073	2233.0	2681.0
	θ[5]	6.608	5.945	-3.666	18.663	0.121	0.085	2269.0	2632.0
	θ[6]	6.774	5.444	-3.065	17.534	0.119	0.084	1959.0	2302.0
	θ[7]	5.028	4.937	-4.818	13.872	0.097	0.07	2386.0	3902.0
	θ[8]	3.742	5.472	-6.317	14.578	0.11	0.078	2294.0	2851.0

As you noticed, the order of θ parameter is quite different between the two.

θ[1] of the chn_post corresponds to θ[5] of the idata.posterior.

Which one is correct?

As arviz package provides this data as zdata=az.load_arviz_data("centered_eight"), this will give us the correct one.

    summarystats(zdata.posterior)

	variable	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail
	mu	4.093	3.372	-2.118	10.403	0.215	0.152	250.0	643.0
	theta[1]	6.026	5.782	-3.707	17.337	0.291	0.206	348.0	743.0
	theta[2]	4.724	4.736	-4.039	13.999	0.199	0.142	471.0	1018.0
	theta[3]	3.576	5.559	-6.779	13.838	0.248	0.175	463.0	674.0
	theta[4]	4.478	4.939	-5.528	13.392	0.199	0.141	503.0	666.0
	theta[5]	3.064	4.642	-5.972	11.547	0.234	0.166	380.0	833.0
	theta[6]	3.821	4.979	-5.507	13.232	0.211	0.15	516.0	1104.0
	theta[7]	6.25	5.436	-3.412	16.92	0.237	0.168	402.0	1026.0
	theta[8]	4.544	5.521	-5.665	15.266	0.231	0.163	449.0	1084.0
	tau	4.089	3.001	0.569	9.386	0.252	0.178	79.0	54.0

This give the conclusion that from_mcmcchains may have disrupted the correct sequence of theta parameters.

What do you think? I'm just confused....

Thanks

[sethaxen]

Thanks for the issue! You're right, that doesn't look good. I just re-ran the Quickstart example on my machine, and summarystats(turing_chns) and summarystats(from_mcmcchains(turing_chns)) both agree for me. Can you give a complete code example that shows this result for me to test with and also post the output of Pkg.status()?

[AnselmJeong]

Thanks for prompt reply...

Here is my code leading to the said inconsistency...
Data were taken from the original arviz homepage.

y = [28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0]
sigma = [15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0]

@model function eight(;y=missing, σ=sigma, N=8)
    y = y===missing ? Array{Float64}(undef, N) : y
    
    μ ~ Normal(0, 5)
    τ ~ truncated(Cauchy(0, 5), 0, Inf)
    θ ~ filldist(Normal(μ, τ), N)
    
    @. y ~ Normal(θ, σ)
end  

model  = eight(;y);
chn_post = sample(model, NUTS(), 1_000);

summarystats(chn_post)

idata = from_mcmcchains(chn_post)
summarystats(idata)

sample() can be a single chain or multiple chains

chn_single_post = sample(model_sample, NUTS(), 1_000);
chn_multiple_post = sample(model_sample, NUTS(), MCMCThreads(), 1_000, 2);

But, this reordering issue applies to both of them..

Pkg.status()

  [cbdf2221] AlgebraOfGraphics v0.4.9
  [131c737c] ArviZ v0.5.5
  [336ed68f] CSV v0.8.5
  [13f3f980] CairoMakie v0.6.3
  [324d7699] CategoricalArrays v0.10.0
  [8be319e6] Chain v0.4.7
  [a93c6f00] DataFrames v1.2.1
  [85a47980] Dictionaries v0.3.10
  [31c24e10] Distributions v0.24.18
  [634d3b9d] DrWatson v2.1.2
  [becb17da] Feather v0.5.9
  [e9467ef8] GLMakie v0.4.4
  [7073ff75] IJulia v1.23.2
  [c7f686f2] MCMCChains v4.13.1
  [ee78f7c6] Makie v0.15.0
  [8b842266] PalmerPenguins v0.1.2
  [91a5bcdd] Plots v1.19.4
  [c3e4b0f8] Pluto v0.15.1
  [7f904dfe] PlutoUI v0.7.9
  [438e738f] PyCall v1.92.3
  [6f49c342] RCall v0.13.12
  [8ce77f84] Soss v0.20.0
  [90137ffa] StaticArrays v1.2.9
  [2913bbd2] StatsBase v0.33.9
  [4c63d2b9] StatsFuns v0.9.8
  [f3b207a7] StatsPlots v0.14.26
  [9e3dc215] TimeSeries v0.23.0
  [fce5fe82] Turing v0.16.6
  [276b4fcb] WGLMakie v0.4.4
  [ade2ca70] Dates

[AnselmJeong]

The weird thing is that, as you inspect the chn_post, the ordering of θ parameters is not as it seems.

group(chn_post, :θ)

Output:

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

Iterations        = 501:1:1500
Number of chains  = 1
Samples per chain = 1000
Wall duration     = 0.5 seconds
Compute duration  = 0.5 seconds
parameters        = θ[5], θ[7], θ[8], θ[2], θ[1], θ[4], θ[6], θ[3]
internals         = 

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

        θ[1]    7.0462    6.4852     0.2051    0.4285   268.2643    1.0079     ⋯
        θ[2]    5.3364    5.1565     0.1631    0.2275   351.0179    0.9996     ⋯
        θ[3]    3.7642    5.6877     0.1799    0.3353   375.4434    0.9991     ⋯
        θ[4]    4.8493    5.3475     0.1691    0.3289   346.8819    1.0046     ⋯
        θ[5]    3.1562    5.1128     0.1617    0.2497   351.7596    1.0020     ⋯
        θ[6]    3.8590    5.2401     0.1657    0.3042   358.8378    1.0002     ⋯
        θ[7]    7.2531    5.9002     0.1866    0.3476   262.3535    1.0127     ⋯
        θ[8]    5.1030    6.0514     0.1914    0.3020   414.4987    1.0040     ⋯

The output indicates that the order of the parameters is θ[5], θ[7], θ[8], θ[2], θ[1], θ[4], θ[6], θ[3] not θ[1], θ[2], θ[3], θ[4], θ[5], θ[6], θ[7], θ[8]. This reordering clearly explains the inconsistency between summarystats(chn_post) and summarystats(from_mcmcchains(chn_post)).

Then the question would be which one is at fault. I want to make sure which θ belongs to which school. The above inspection suggests that summarystats(chn_post) are not reporting properly.

Do you have any idea?

Thanks again...

[AnselmJeong]

I think I found the possible reason.

Comment from the maintainer of MCMCChains

I raised this issue in MCMCChains github page.

But I'm still not sure which one is confusing the order.

[sethaxen]

Yes, that's absolutely right. Our code assumes that the entries of the variables are sorted in a natural order, but it looks like with Turing v0.16, that doesn't seem to be the case anymore. I'll push a fix so we no longer make assumptions about sorting.