Plotting credible intervals in Julia from Turing model

Question

Ok so I figured out how to plot the credible intervals for a univariate linear model in Turing.jl using the following code (I'm replicating Statistical rethinking by McElreath) This particular exercise is in chapter 4. If anyone has already plotted these types of models with Turing and can give me a guide, it would be great!!!

Univariate model code:

    using Turing
    using StatsPlots
    using Plots

    height = df2.height

    weight = df2.weight

    @model heightmodel(y, x) = begin
        #priors
        α ~ Normal(178, 100)
        σ ~ Uniform(0, 50)
        β ~ LogNormal(0, 10)

        x_bar = mean(x)
        #model

            μ = α .+ (x.-x_bar).*β
        y ~ MvNormal(μ, σ)
    end

    chns = sample(heightmodel(height, weight), NUTS(), 100000)

    ## code 4.43
    describe(chns) |> display

    # covariance and correlation

    alph = get(chns, :α)[1].data

    bet = get(chns, :β)[1].data

    sigm = get(chns, :σ)[1].data


    vecs = (alph[1:352], bet[1:352])
    arr = vcat(transpose.(vecs)...)'

    ss = [vec(alph + bet.*(x)) for x in 25:1:70]

    arrr = vcat(transpose.(ss)...)'

    plot([mean(arrr[:,x]) for x in 1:46],25:1:70, ribbon = ([-1*(quantile(arrr[:,x],[0.1,0.9])[1] - mean(arrr[:,x])) for x in 1:46], [quantile(arrr[:,x],[0.1,0.9])[2] - mean(arrr[:,x]) for x in 1:46]))

Credible interval Univariate:

However, when I try to replicate it with a multivatiate function, very strange things are drawn:

Multivariate model code:

    weight_s = (df.weight .-mean(df.weight))./std(df.weight)

    weight_s² = weight_s.^2

    @model heightmodel(height, weight, weight²) = begin
        #priors
        α ~ Normal(178, 20)
        σ ~ Uniform(0, 50)
        β1 ~ LogNormal(0, 1)
        β2 ~ Normal(0, 1)
        #model
        μ = α .+ weight.*β1 + weight².*β2
        height ~ MvNormal(μ, σ)
    end

    chns = sample(heightmodel(height, weight_s, weight_s²), NUTS(), 100000)

    describe(chns) |> display

    ### painting the fit

    alph = get(chns, :α)[1].data

    bet1 = get(chns, :β1)[1].data

    bet2 = get(chns, :β2)[1].data



    vecs = (alph[1:99000], bet1[1:99000], bet2[1:99000])
    arr = vcat(transpose.(vecs)...)'


    polinomial = [vec(alph + bet1.*(x) + bet2.*(x.^2)) for x in -2:0.01:2]

    arrr = vcat(transpose.(polinomial)...)'

    plot([mean(arrr[:,x]) for x in 1:401],-2:0.01:2, ribbon = ([-1*(quantile(arrr[:,x],[0.1,0.9])[1] - mean(arrr[:,x])) for x in 1:46], [quantile(arrr[:,x],[0.1,0.9])[2] - mean(arrr[:,x]) for x in 1:46]))

Credible interval Univariate:

user_8683182 · Answer

In the Julia Slack channel (https://slackinvite.julialang.org/) Jens was kind enough to give me the answer. The credit goes to him --He doesn't have a SO account.

The main problem was that I was lacking simplicity and was trying to plot the mean the wrong way --in a very inefficient and strange way might I add. Each parameter has a vector which corresponds to 99000 draws from the posterior distribution, I was trying to draw the mean through the matrix but it's much easier to define the median first and then plot it, then you don't make mistakes as I did calculating the mean.

old code

vecs = (alph[1:99000], bet1[1:99000], bet2[1:99000])
arr = vcat(transpose.(vecs)...)'
[mean(arrr[:,x]) for x in 1:401]

can be written as:

testweights = -2:0.01:2
arr = [fheight.(w, res.α, res.β1, res.β2) for w in testweights]
m = [mean(v) for v in arr]

Moreover, the way Jens defined the Credible intervals is much more elegant and Julianic:

Jens' code:

quantiles = [quantile(v, [0.1, 0.9]) for v in arr]
lower = [q[1] - m for (q, m) in zip(quantiles, m)]
upper = [q[2] - m for (q, m) in zip(quantiles, m)]

My code:

ribbon = ([-1*(quantile(arrr[:,x],[0.1,0.9])[1] - mean(arrr[:,x])) for x in 1:46], [quantile(arrr[:,x],[0.1,0.9])[2] - mean(arrr[:,x]) for x in 1:46]))

Complete Solution:

weight_s = (d.weight .-mean(d.weight))./std(d.weight)
height = d.height
@model heightmodel(height, weight) = begin
    #priors
    α ~ Normal(178, 20)´
    σ ~ Uniform(0, 50)
    β1 ~ LogNormal(0, 1)
    β2 ~ Normal(0, 1)
    #model
    μ = α .+ weight .* β1 + weight.^2 .* β2
    # or μ = fheight.(weight, α, β1, β2) if we are defining fheigth anyways
    height ~ MvNormal(μ, σ)
end
chns = sample(heightmodel(height, weight_s), NUTS(), 10000)
describe(chns) |> display
res = DataFrame(chns)
fheight(weight, α, β1, β2) = α + weight * β1 + weight^2 * β2
testweights = -2:0.01:2
arr = [fheight.(w, res.α, res.β1, res.β2) for w in testweights]
m = [mean(v) for v in arr]
quantiles = [quantile(v, [0.1, 0.9]) for v in arr]
lower = [q[1] - m for (q, m) in zip(quantiles, m)]
upper = [q[2] - m for (q, m) in zip(quantiles, m)]
plot(testweights, m, ribbon = [lower, upper])

Plotting credible intervals in Julia from Turing model

1 Answers