# Parameter estimation problem in Turing
**URL:** https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147
**Category:** Probabilistic Programming
**Tags:** turing
**Created:** [July 13, 2022, 6:33am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147 "2022-07-13T06:33:20Z")
**Posts on this page:** 12
**Page:** 1
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 13, 2022, 6:33am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/1 "2022-07-13T06:33:20Z")
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.
```julia
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:
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!
---
### Author:  [@wc4wc4wc4](https://discourse.julialang.org/u/wc4wc4wc4)
#### Post date: [July 13, 2022, 8:01am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/2 "2022-07-13T08:01:27Z")
Hi again,
Try to make a [plot of the chains](https://beta.turing.ml/MCMCChains.jl/dev/statsplots/), that’s usually very useful for understanding your results.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 13, 2022, 8:13am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/3 "2022-07-13T08:13:53Z")
Thank you for the suggestion!
I will add the visualization for the chain when doing MCMC.
---
### Author:  [@skleinbo](https://discourse.julialang.org/u/skleinbo)
#### Post date: [July 13, 2022, 8:33am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/4 "2022-07-13T08:33:32Z")
At first glance, the model you have written is different from the generative process. In the latter, each observation y\_i has its own \beta\_i associated. While in the model, one \beta is sampled which informs all y .
Instead of one \beta, I believe you want something like
```julia
β ~ filldist(Gamma(α0, β0), N)
for i in 1:N
y[i] ~ Gamma(α,β[i])
end
```
which for 1000 samples along one chain yields for me
```julia
Summary Statistics
parameters mean std naive_se mcse ess rhat ess_per_sec
Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64
α 0.4783 0.1643 0.0052 0.0232 26.4629 1.0126 0.9557
α0 0.5508 0.2129 0.0067 0.0302 25.4590 1.0162 0.9194
β0 0.0467 0.0191 0.0006 0.0016 121.8532 1.0064 4.4006
```
That’s not _too_ bad I’d say.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 13, 2022, 8:39am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/5 "2022-07-13T08:39:11Z")
@skleinbo Thank you very much for providing this solution.
Yes, your implementation of the model is exactly what I want.
I will further develop based on your solution.
---
### Author:  [@skleinbo](https://discourse.julialang.org/u/skleinbo)
#### Post date: [July 13, 2022, 8:46am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/6 "2022-07-13T08:46:44Z")
You may probably also want to widen the priors (which I did!). For example, the prior for \alpha\sim\text{Gamma}(0.1,0.2) has P(0.3\<\alpha\<0.4) \approx 0.005.
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 13, 2022, 8:47am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/7 "2022-07-13T08:47:28Z")
Hi, again. Can you help to provide the whole example of your implementation?
I just run into some error saying:
```julia
Output exceeds the [size limit](command:workbench.action.openSettings?[). Open the full output data [in a text editor](command:workbench.action.openLargeOutput?74549bc6-1cc2-4911-a2b2-89cd237c0afe)
BoundsError: attempt to access 0-element Vector{Float64} at index [1]
```
Thank you!
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 13, 2022, 8:48am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/8 "2022-07-13T08:48:42Z")
Yes. I will try this suggestion.
Thank you very much!
---
### Author:  [@skleinbo](https://discourse.julialang.org/u/skleinbo)
#### Post date: [July 13, 2022, 8:50am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/9 "2022-07-13T08:50:36Z")
Note that I have changed the model to accept y directly, instead of conditioning it later like you did.
```julia
julia> using Turing, MCMCChains, Distributions
using StatsPlots
# construct the probabilistic model
@model function rdeg2(y)
N = length(y)
# define prior
α ~ Gamma(1.0, 0.2)
α0 ~ Gamma(1.0, 0.2)
β0 ~ Gamma(1.0, 0.2)
# Likelihood
β ~ filldist(Gamma(α0, β0), N)
for i in 1:N
y[i] ~ Gamma(α,β[i])
end
return y
end
# collect data
data = Float64[]
for i in 1:120
β = rand(Gamma(0.6, 0.04))
push!(data, rand(Gamma(0.4, β)))
end
model = rdeg2(data)
chain = sample(model, NUTS(), 1000)
```
It doesn’t make a difference for the inference, but if you define `rdeg(;N)` like you originally did, you’ll need to define
```julia
y = Vector{Float64}(undef, N) # needs to be declared if not argument of the model
for i in 1:N
y[i] ~ Gamma(α,β[i])
end
```
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 13, 2022, 9:48am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/10 "2022-07-13T09:48:19Z")
@skleinbo Thank you very much for the solution and explanation.
There is one thing I need to ask. I run the program and got
Here the parameters include all βs?
---
### Author:  [@skleinbo](https://discourse.julialang.org/u/skleinbo)
#### Post date: [July 14, 2022, 8:15am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/11 "2022-07-14T08:15:58Z")
Yes, because as I said, each observation comes from a distribution with its own \beta\_i.
Had you drawn observations like
```julia
\beta = rand(Gamma(\alpha0,\beta0)
y = rand(Gamma(\alpha,\beta), 120)
```
your original model would be the right one.
Just for fun, you could plot a histogram of all the \beta samples from the chain(s) and compare to the generating distribution \text{Gamma}(0.6,0.04).
---
### Author:  [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)
#### Post date: [July 14, 2022, 10:16am UTC](https://discourse.julialang.org/t/parameter-estimation-problem-in-turing/84147/12 "2022-07-14T10:16:35Z")
Thank you very much for the reply!