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."

[torfjelde]

I'm very excited about this 👀

I have a few minor comments, but will approve as it looks good:)

[devmotion]

It seemed you wanted to push a few additional updates @sethaxen? Or is this PR ready from your side?

[sethaxen]

Yeah it's pending a breaking changes to PosteriorStats I need to finish up, to avoid another breaking release here.

[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.