StatisticalRethinkingJulia / StatisticalRethinking.jl

Julia package with selected functions in the R package `rethinking`. Used in the SR2... projects.
MIT License
386 stars 32 forks source link

Plotbounds for Turing #94

Closed orangebacked closed 4 years ago

orangebacked commented 4 years ago

Is there any code to plot the bounds using Turing? I have succeeded in plotting the linear models but haven't been able to plot the quadratic or cubic models.

goedman commented 4 years ago

You mean using the plotbounds() method? Did you try converting the MCMCChains.Chains object to a DataFrame, e.g. DataFrame(chns, append_chains=true)?

If not, can you send me a MWE and I'll take a look. Clearly we're just starting creating a Turing version of StatisticalRethinking.jl.

orangebacked commented 4 years ago

Dear Mr. Goedman,

I was referring indeed to the plotbounds() method. I managed to plot it with the help of the Julia of the Julia community in Slack (you can find the discussion here https://stackoverflow.com/questions/62028147/plotting-credible-intervals-in-julia-from-turing-model/62070008?noredirect=1#comment109792812_62070008). This is the code (I hope it is useful):

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

I'm transcribing the book with julia and Turing.jl, albeit slowly; let me know if I can colaborate! I saw this message today (I'm sorry for the delay), I guess I have to activate the notifications in my email.

Kind Regards, Emiliano Isaza V

goedman commented 4 years ago

Emiliano,

I apologize for not responding to your post. I don't understand how I missed your post 10 days ago, but somehow I must have deleted the notification.

I'll go over above code today and see where it differs from plotbounds() in SR and if it shows ways to improve plotbounds().

Best, Rob

goedman commented 4 years ago

Yes, very nice. I have started to use ribbons in other contexts as well but should probably rework this into plotbounds(). The other part I would like to allow in plotbounds() is to pass in a user defined link function like fweight() in your code above.

Have you looked at @karajan9 's work on TuringModels and SR? Looks really neat!