Построение вероятных интервалов в Джулии из модели Тьюринга

Question

Хорошо, поэтому я понял, как построить вероятные интервалы для одномерной линейной модели в Turing.jl, используя следующий код (я копирую «Статистическое переосмысление Макэлрита»). Это конкретное упражнение находится в главе 4. Если кто-то уже построил график эти типы моделей с Тьюрингом и могут дать мне руководство, было бы здорово !!!

Одномерный код модели:

    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]))

Достоверный интервал Одномерный:

Однако, когда я пытаюсь воспроизвести его с помощью многомерной функции, рисуются очень странные вещи:

Код многомерной модели:

    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]))

Достоверный интервал Одномерный:

Answer 1

В канале Julia Slack (https://slackinvite.julialang.org/) Йенс был достаточно любезен, чтобы дать мне ответ. Заслуга принадлежит ему - у него нет учетной записи SO.

Основная проблема заключалась в том, что мне не хватало простоты, и я пытался неправильно изобразить среднее - очень неэффективным и странным способом, могу добавить. У каждого параметра есть вектор, который соответствует 99000 dr aws из апостериорного распределения, я пытался провести среднее значение через матрицу, но намного проще сначала определить медиану, а затем построить ее, тогда вы не ошибетесь, поскольку Я вычислил среднее значение.

старый код

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

можно записать как:

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

Более того, способ, которым Йенс определил достоверные интервалы, намного элегантнее и Джулиани c:

Код Йенса:

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)]

Мой код:

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]))

Полное решение:

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])