Using Turing for MCMC estimation

General Usage Performance
,
1

Hi,
I am new to Turing PPL.
I have a random variable X which is Gamma distributed, i.e., X ~Ga(α, β). And the scale β parameter of X is also obeys a Gamma law, namely, β ~ Ga(α0, β0). Then I want to infer the above parameters from the observed data of X. I was thinking the parameters are α, α0, β0.
I implemented the following solution based on the NUTS sampler in Turing.

using Turing, MCMCChains, Distributions
using StatsPlots
# construct the probabilistic model
@model function rdeg(; N::Int)
    # define prior
    α ~ Gamma(0.1, 0.2)
    α0 ~ Gamma(0.1, 0.2)
    β0 ~ Gamma(0.1, 0.2)

    # Likelihood
    β ~ Gamma(α0, β0)
    y ~ filldist(Gamma(α, β), N)
    # y ~ filldist(Gamma(α, Gamma(α0, β0)), N)
    return y

end
rdeg(y::AbstractVector{<:Real}) = rdeg(; N=length(y)) | (; y);

# collect data
data = Float64[]
for i in 1:120
    β = rand(Gamma(0.6, 0.04))
    push!(data, rand(Gamma(0.4, β)))
end
model = rdeg(data)
# sampler = HMC(0.09, 15)
sampler = NUTS() # a nice sampler
chain = sample(model, sampler, MCMCSerial(), 10_000, 4, progress=false)

As seen in above program, I have set α, α0, β0 = 0.4, 0.6, 0.04; And I got output of:

turing
turing883×764 96 KB

First, I see the estimation of α, α0, β0 are not so accurate, α is nearly Ok, but α0, β0 are not.
The second problem is I am not sure if my implementation of inference model in @model block is correct.
I see in the results, the parameters are four: α, β, α0, β0.

Thank you very much for checking my problems!

2

This program is not working, how can obtain the estimations for all parameters?

You’re more likely to get a helpful reply if you provide more details. What does “not working” mean? Is there an error message? If yes, provide the message. Is there no error, but the output is incorrect? What is the output? How do you know it’s incorrect?

Without that, I’ll guess at one problem: data = [] creates a Vector{Any}, but you’ve define rdeg to only work with AbstractVector{<:Real}. Initialize data = Float64[] instead should fix that problem.

2 Likes
3

Thank you very much Paul!
Your reply solves the error in my program, now I can obtain the
results.

4

Also, for future help with Turing, the performance section is not really the best place to ask. You would probably get more answers by asking here: Probabilistic programming - JuliaLang.

1 Like
5

@wc4wc4wc4 Thank you very much for sharing this information!
I will try on the link you shared.

1 Like