Import and extend PosteriorStats

[sethaxen]

This PR makes the following replacements everywhere:

• hpd -> PosteriorStats.hdi (hpd retains its old behavior and is now deprecated)
• summarize -> PosteriorStats.summarize
• ChainDataFrame -> PosteriorStats.SummaryStats

The replacement of ChainDataFrame and the slight change in API and behavior of the methods makes this a breaking change.

Implements #430

e.g.

julia> val = rand(500, 2, 3);

julia> chn
Chains MCMC chain (500×2×3 Array{Float64, 3}):

Iterations        = 1:1:500
Number of chains  = 3
Samples per chain = 500
parameters        = param_1, param_2


Summary Statistics
          mean    std    hdi_3%  hdi_97%  mcse_mean  mcse_std  ess_tail  ess_bulk  rhat 
 param_1  0.50  0.283  0.00745     0.929     0.0073    0.0034      1480      1531  1.00
 param_2  0.49  0.285  0.000140    0.926     0.0074    0.0034      1499      1506  1.00

Quantiles
            2.5%  25.0%  50.0%  75.0%  97.5% 
 param_1  0.0252  0.254  0.510  0.745  0.969
 param_2  0.0180  0.242  0.483  0.727  0.969

julia> hdi(chn)
HDI
             lower  upper 
 param_1  0.00745   0.929
 param_2  0.000140  0.926

julia> describe(chn)
2-element Vector{SummaryStats}:
 Summary Statistics (2 rows, 10 cols)
 Quantiles (2 rows, 6 cols)

[sethaxen]

Lacking a parameter column meant the show method was broken. But since there is a changerate for every parameter, it makes more sense to do the same thing as gelmandiag_multivariate and return a SummaryStats for the marginal values and return the multivariate changerate separately.

[torfjelde]

Genuine question: what is the change-rate in this context?

[sethaxen]

From inspecting the code, it's, for each parameter and chain, the fraction of draws that are different from the previous draw. I suppose it's similar to "acceptance rate."

[sethaxen]

@devmotion @torfjelde This is ready for final review.

[sethaxen]

Currently PosteriorStats has no try/catch mechanism for if a given statistic fails. That causes chains with less than 10 draws to error upon computation of summary stats (or display). For PosteriorStats that makes sense, but it's a major nuisance for MCMCChains. I actually think this has come up enough that it makes sense for MCMCDiagnosticTools to raise an informative warning and return NaNs in such cases. Will open an issue there.

[devmotion]

Oh, I missed #431 (comment), probably due to the Christmas break. I added the PR to my todo list 🙂

[devmotion]

This seems a bit inefficient (particular the use of stack(map(stack) - maybe we should preallocate the desired array and copy directly? Possibly a utility function could be added to PosteriorStats?

[sethaxen]

It may not be the most efficient, but compared to the actual autocor call, it should not be the computational bottleneck, and not using stack requires substantially more code.

It might make sense to add to PosteriorStats a function that converts a SummaryStats to a matrix, but this would effectively do the same thing as the inner stack. What kind of utility function do you have in mind?

[devmotion]

What kind of utility function do you have in mind?

Basically a convert, Array/Matrix constructor, or other function that performs the same thing (if the first two alternatives are not desirable). So that downstream packages (such as MCMCChains) don't have to deal with internals of SummaryStats.

[sethaxen]

This line is not using any internals of SummaryStats, whose docstring explains it is an OrderedDict-like Table that can be iterated/indexed by its column names.

I think the methods we're missing in PosteriorStats are

function Base.getindex(stats::SummaryStats, cols::Union{Colon,AbstractVector{Int},AbstractVector{Symbol}})
    cols isa Colon && return Tables.matrix(stats)
    return stack(Tables.getcolumn(stats, k) for k in cols)
end
Base.firstindex(s::SummaryStats) = 1
Base.lastindex(s::SummaryStats) = length(s)

Then this would be

stack(ac_i[2:end] for ac_i in ac)

but this is no more efficient.

[devmotion]

This line is not using any internals of SummaryStats, whose docstring explains it is an OrderedDict-like Table that can be iterated/indexed by its column names.

I was referring to the 2:end which seems a bit special. I missed that the docstring mentions that the first column is reserved for parameter names.

[devmotion]

Maybe a good opportunity to make the code a bit more efficient:

Suggested change

	chain_dic = Dict(zip(quantile(chains)[2], quantile(chains)[5]))
	quantile_chains = quantile(chains)
	chain_dic = Dict(zip(quantile_chains[2], quantile_chains[5]))

(it would also be nice if it would be possible to use something more descriptive than 2 and 5)

[sethaxen]

Hmm, yeah this and the code immediately below are a bit of a mess. Will rework to use the new methods.

[sethaxen]

Yeah it was able to be simplified a lot. But also, for n parameters, it was computing the median n^2 times. This has now been fixed.

[torfjelde]

Added a few comments!

Am very excited about this though; bloody great stuff @sethaxen ❤️

[torfjelde]

What's the meaning of "unimodal" here?

[sethaxen]

The default estimator used by hdi assumes that the distribution is unimodal. The alternative (not yet in PosteriorStats but doable with HighestDensityRegions.jl) for multimodal distributions first fits a KDE and then partitions into one or more intervals by density. For fast summary statistics, the unimodal version is more useful.

[torfjelde]

Genuine question: what is the change-rate in this context?

[sethaxen]

This is ready for another review. Also, it would be nice to get input on arviz-devs/PosteriorStats.jl#25, since that's relevant for this PR.

[yebai]

@penelopeysm can you help review this PR, and get it merged?

[devmotion]

I thought arviz-devs/PosteriorStats.jl#25 holds back this PR?

[sethaxen]

@devmotion is right, I'd like to merge that breaking change and also some breaking changes to make the SummaryStats object easier to interact with (and a little more like ChainDataFrame). Otherwise when I make those changes to PosteriorStats, it would require another major version bump here and a minor one for Turing.

After a long break I'm back to working on arviz PRs so I should be able to wrap that up in the next couple of weeks.