# Working with posteriors

**URL:** <https://discourse.julialang.org/t/working-with-posteriors/19846>\
**Category:** Statistics\
**Tags:** question, splitapplycombine\
**Created:** [January 20, 2019, 6:13pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846 "2019-01-20T18:13:15Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 20, 2019, 6:13pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/1 "2019-01-20T18:13:15Z")

</div>

Finding that I am writing a lot of repetitive code, I am trying to rethink the way I work with posterior distributions (eg from MCMC).

I can typically reduce most operations on posteriors to _transformations_, which operate on a particular draw from the posterior, and _elementwise_ operations, eg calculating means and quantiles for various, potentially nested containers, eg an upper-triangular matrix

For transformations, the best container type is an `AbstractVector{T}` or similar, where `T` can be a `struct` or a `NamedTuple` of various scalars, arrays, etc; transforming to something similar (or a `Tuple`) one draw at a time. This is easy to program and iterate. Typical examples would be posterior predictive checks, transformed parameters, and summary statistics (eg RMSE and equivalents).

On the other hand, elementwise operations typically work best on something `AbstractMatrix`-like, where one can map columns using `eachcol` or `mapslices`.

One typically wants to convert from one representation to another multiple times.

I am curious how other Julia users approach this, and what existing libraries would be recommended. Packages with somewhat similar (but not identical) goals that I am aware of include

> **[GitHub - SciML/RecursiveArrayTools.jl: Tools for easily handling objects like...](https://github.com/SciML/RecursiveArrayTools.jl)**
>
> Tools for easily handling objects like arrays of arrays and deeper nestings in scientific machine learning (SciML) and other applications - GitHub - SciML/RecursiveArrayTools.jl: Tools for easily h...

> **[GitHub - JuliaArrays/StructArrays.jl: Efficient implementation of struct...](https://github.com/JuliaArrays/StructArrays.jl)**
>
> Efficient implementation of struct arrays in Julia - GitHub - JuliaArrays/StructArrays.jl: Efficient implementation of struct arrays in Julia

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [January 20, 2019, 6:23pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/2 "2019-01-20T18:23:53Z")

</div>

Maybe [Queryverse](http://www.david-anthoff.com/jl4ds/stable/)? As I understand it abstracts this concept and aims to easily connect a wide variety of data sources and sinks

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 20, 2019, 6:37pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/3 "2019-01-20T18:37:20Z")

</div>

Maybe I was not clear (my bad). To make things concrete, consider the mock example

```julia
posterior_draws = [(a = i, b = fill(i, 3), c = fill(i, 3, 3)) for i in 1:10]
_flatten(d) = vcat(d.a, d.b, vec(d.c))
posterior_matrix = collect(permutedims(mapreduce(_flatten, hcat, posterior_draws)))

```

Sometimes you want to work with draws, eg

```julia
function do_stuff(d) # operate on posterior draw
    (statistic = d.a + sum(d.c), )
end
map(do_stuff, posterior_draws)

```

Then sometimes you prefer elementwise:

```julia
map(mean, eachcol(posterior_matrix))

```

My goal is to go back and forth between the two representations in some organized manner. This probably involves defining some kind of schema for the matrix representation’s column layout. I can do this, but wanted to ask before making yet another package.

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [January 20, 2019, 6:59pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/4 "2019-01-20T18:59:09Z")

</div>

Is the current problem you’re hitting mostly trouble with performance, or ease of use? Where do you find the packages you linked to missing the mark?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 20, 2019, 7:09pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/5 "2019-01-20T19:09:48Z")

</div>

I am going for ease of use. Mosly, I have \<10k posterior draws, so programming them is the costly part, not the runtime.

The packages I linked do something _similar conceptually_, but not what I describe above (or maybe I am using them wrong — in which case rewriting the examples above would be helpful).

---

<div class="post-metadata">

**Author:** ![cscherrer](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/cscherrer/32/7631_2.png) [@cscherrer](https://discourse.julialang.org/u/cscherrer)\
**Post date:** [January 20, 2019, 7:49pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/6 "2019-01-20T19:49:06Z")

</div>

I don’t know that this is the best approach, but I’ve usually tried to keep things as an iterable over observations, and then do things like `getfield.(posterior_draw, :a)` when I need a column.

There seem to be quite a few packages trying to do what you describe or very similar things. [DataStreams.jl](https://github.com/JuliaData/DataStreams.jl) is another worth checking out, as is [IndexedTables.jl](https://github.com/JuliaComputing/IndexedTables.jl) that’s used by JuliaDB.

Seems like one potential difficulty is that your columns have some structure - hierarchical, even if only one level of hierarchy.

---

<div class="post-metadata">

**Author:** ![piever](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/piever/32/1815_2.png) [@piever](https://discourse.julialang.org/u/piever)\
**Post date:** [January 20, 2019, 8:51pm UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/7 "2019-01-20T20:51:38Z")

</div>

IndexedTables is “powered” by StructArrays (an `IndexedTable` is mainly a wrapped `StructArray` with information on which columns are sorted). Other than the packages you mention (StructArrays and RecursiveArrayTools), I think there is a simple idea that may be worth exploring:  
StructArrays accepts arbitrary abstract arrays as columns, so you can combine it with a custom array types that is stored as an N-dimensional array but is an `AbstractVector`:

```julia
using StructArrays
using StructArrays: fieldarrays

struct FakeVector{T, N, M} <: AbstractVector{Array{T, M}}
    array::Array{T, N}
    function FakeVector(array::Array{T, N}) where {T, N}
        new{T, N, N-1}(array)
    end
end

eltypelength(f::FakeVector{T, N, M}) where {T, N, M} = M
Base.size(f::FakeVector) = (size(f.array, 1),)
Base.getindex(f::FakeVector, i) = f.array[i, (Colon() for x in 1:eltypelength(f))...] 
function Base.setindex!(f::FakeVector, val, i)
    f.array[i, (Colon() for _ in 1:eltypelength(f))...] = val
end
reshape2D(v::AbstractArray) = reshape(v, size(v, 1), :)
reshape2D(f::FakeVector) = reshape2D(f.array)

# usage
N = 1000
vals = rand(N)
gradients = FakeVector(rand(N, 3))
hessians = FakeVector(rand(N, 3, 3))
draws = StructArray(vals = vals, gradients = gradients, hessians = hessians)
draws_matrix = mapreduce(reshape2D, hcat, fieldarrays(draws))

# accessing individual arrays:
draws.gradients.array

```

In short I think that something like the `FakeVector` type above + StructArrays would solve it and I actually would like a `FakeVector` type like that. I may go as far as putting it in StructArrays or IndexedTables as it really helps by allowing to use tabular data machinery when working with arrays of more than 1 dimension. Does a package implementing it exist already?

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [January 22, 2019, 10:45am UTC](https://discourse.julialang.org/t/working-with-posteriors/19846/8 "2019-01-22T10:45:13Z")

</div>

After a bit of experimentation, I found that I can mostly work with vectors of things (`NamedTuple`s, matrices, etc) and `map`, occasionally manipulating the arrays with @andyferris’s excellent

> **[GitHub - JuliaData/SplitApplyCombine.jl: Split-apply-combine strategies for...](https://github.com/JuliaData/SplitApplyCombine.jl)**
>
> Split-apply-combine strategies for Julia. Contribute to JuliaData/SplitApplyCombine.jl development by creating an account on GitHub.

and using @mauro3’s

> **[GitHub - mauro3/Parameters.jl: Types with default field values, keyword...](https://github.com/mauro3/Parameters.jl)**
>
> Types with default field values, keyword constructors and (un-)pack macros - GitHub - mauro3/Parameters.jl: Types with default field values, keyword constructors and (un-)pack macros

to make working with `NamedTuple` very convenient. A worked example is below.

```julia
using SplitApplyCombine, Parameters, Statistics, Distributions, PGFPlotsX

###
### generate data (with fatter tails) and toy example posterior.
###
### (The posterior does not come from any valid inference, to make
### the example self-contained. It is not valid statistical methodology).

data = rand(TDist(3), 100); # fat tails, will be detected by PPC
posterior = [(μ = randn()/10, σ = rand(Uniform(1.5, 2))) for _ in 1:1000];

# simulate data from posterior
sim = let N = length(data)
    map(posterior) do p
        @unpack μ, σ = p
        rand(Normal(μ, σ), N)
    end
end

function stats(d, p)
    @unpack μ, σ = p
    q5, q25, q75, q95 = quantile(@.((d - μ)/σ), [0.05, 0.25, 0.75, 0.95])
    [q75 - q25, q95 - q5]
end

data_stats = combinedims(map(p -> stats(data, p), posterior))
sim_stats = combinedims(map(stats, sim, posterior))

function makeplot(varlabel, data_stat, sim_stat)
    ab = [extrema(vcat(data_stat, sim_stat))...]
    @pgf Axis({ xlabel = "data " * varlabel, ylabel = "simulated " * varlabel },
              Plot({ blue }, Table(ab, ab)), # 45° line
              Plot({ only_marks, mark = "o", opacity = 0.5 }, Table(data_stat, sim_stat)))
end

makeplot(raw"standardized $q_{75} - q_{25}$", data_stats[1, :], sim_stats[1, :])
makeplot(raw"standardized $q_{95} - q_{5}$", data_stats[2, :], sim_stats[2, :])

```

 ![1](https://global.discourse-cdn.com/julialang/original/3X/9/2/92f6253e2d8d2a82f4c07219ec3679097aaff3a7.png)

 ![2](https://global.discourse-cdn.com/julialang/original/3X/d/5/d545971d26d201bce23351d5cdda971b45ba80c9.png)

These very modular tools make PPC very convenient in Julia. The only missing piece is getting the `NamedTuple`s from the flat matrix data, for which I am working on an experimental package.